I. PLOTS OF TAXONOMIC COMPOSITION

Loading libraries

library(phyloseq)
library(ggplot2)
library(vegan)
library(picante)
library(edgeR)
library("RColorBrewer")
library(scales)
library(grid)
library(reshape2)
library(scales)
library(viridis)
library(hrbrthemes)
library(tidyverse)
library(VennDiagram)
#Installing edgeR

if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")

BiocManager::install("edgeR")

PLOT RELATIVE ABUNDANCE AT THE PHYLUM LEVEL

Load data

metadata <- read.csv(file = "../data/metadata.csv", header = TRUE, row.names = 1)
asv_table <- read.csv("../data/feature_table.csv", header = TRUE, row.names = 1)
taxonomy <- read.csv("../data/taxonomy.csv", header = TRUE, row.names = 1)

Create phyloseq objects

SAM <- sample_data(metadata)
TAX <- tax_table(as.matrix(taxonomy))
ASV <- otu_table(asv_table, taxa_are_rows=TRUE)  
physeq <- merge_phyloseq(ASV, TAX, SAM)

Data visualizations

sample_names(physeq)
##  [1] "Sg.F.D.1"  "Sg.F.D.2"  "Sg.M.D.3"  "Sg.M.D.4"  "Sg.F.D.5"  "Sg.F.D.6" 
##  [7] "Sg.F.D.7"  "Sg.F.D.8"  "Sg.F.D.9"  "Sg.M.D.10" "Sg.F.D.11" "Sg.M.D.12"
## [13] "Sg.M.D.13" "Sg.M.D.14" "Sg.M.D.15" "Ss.M.D.16" "Ss.M.D.17" "Ss.M.D.18"
## [19] "Ss.F.D.19" "Ss.M.D.20" "Ss.F.D.21" "Ss.M.D.22" "Ss.M.D.23" "Ss.M.D.24"
## [25] "Ss.F.D.25" "Ss.M.D.26" "Sa.M.D.27" "Sa.F.D.28" "Sa.M.D.29" "Sa.F.D.30"
## [31] "Sa.M.D.31" "Sa.M.D.32" "Sa.F.D.33" "Sa.F.D.35" "Sa.M.D.36" "Sb.F.D.37"
## [37] "Sb.M.D.38" "Sb.F.D.39" "Sb.M.D.41" "Sb.M.D.42" "Sb.M.D.43" "Sb.F.D.44"
## [43] "Sb.F.D.45" "Sb.M.D.46" "Ss.M.R.52" "Ss.F.R.53" "Ss.M.R.54" "Ss.M.R.55"
## [49] "Ss.F.R.56" "Ss.F.R.57" "Ss.F.R.58" "Ss.M.R.59" "Ss.M.R.60" "Ss.M.R.61"
## [55] "Sg.M.R.62" "Sg.M.R.63" "Sg.M.R.64" "Sg.M.R.65" "Sg.M.R.66" "Sg.F.R.67"
## [61] "Sg.F.R.68" "Sg.F.R.69" "Sg.M.R.70" "Sg.M.R.71" "Sg.F.R.72" "Sg.M.R.74"
## [67] "Sg.M.R.75" "Sa.F.R.72" "Sa.F.R.73" "Sa.M.R.74" "Sa.M.R.75" "Sa.F.R.76"
## [73] "Sa.M.R.77" "Sa.M.R.78" "Sa.M.R.79" "Sa.F.R.80" "Sa.M.R.81" "Sa.F.R.82"
## [79] "Sa.F.R.83" "Sa.F.R.84" "Sb.M.R.85" "Sb.M.R.86" "Sb.M.R.87" "Sb.F.R.88"
## [85] "Sb.F.R.89" "Sb.F.R.90" "Sb.M.R.91" "Sb.F.R.92" "Sb.M.R.93" "Sb.M.R.94"
## [91] "Sb.M.R.95" "Sb.F.R.96" "Sb.F.R.97" "Sb.M.R.98" "Sb.F.R.99"
rank_names(physeq)
## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"
sample_variables(physeq)
##  [1] "Season"               "Species_Scel"         "Sex"                 
##  [4] "SeqDepth"             "Weight"               "SVL"                 
##  [7] "TL"                   "Tb"                   "Ta"                  
## [10] "Ts"                   "Elevation"            "Location"            
## [13] "Date"                 "BarcodeSequence"      "LinkerPrimerSequence"
# Convert to relative abundances 
relative <- transform_sample_counts(physeq = physeq, 
                                    function(ASV) ASV / sum(ASV))

RELATIVE ABUNDACE BY TAXONOMIC LEVEL (MEAN AND STANDARD DEVIATION)

Calculations

data <-
  physeq %>%
  tax_glom("Genus") %>%
  transform_sample_counts(function(x)100* x / sum(x)) %>%
  psmelt() %>%
  as_tibble()

PLOT RELATIVE ABUNDANCE AT THE PHYLUM LEVEL

Create color palette

paleta <- c(brewer.pal(12, "Paired")[1:12], "gray", "orange")
print(paleta)
##  [1] "#A6CEE3" "#1F78B4" "#B2DF8A" "#33A02C" "#FB9A99" "#E31A1C" "#FDBF6F"
##  [8] "#FF7F00" "#CAB2D6" "#6A3D9A" "#FFFF99" "#B15928" "gray"    "orange"

Species + Season

Rel_Abun_Phyl_Season <- plot_bar(physeq = relative, "Sample", fill = "Phylum") + 
  facet_grid(~Species_Scel + Season, space = "free", scales = "free") +
  labs(y="Relative abundance") +
  geom_bar(stat = "identity", pisition = "stack", res = 300) +
  scale_fill_manual(values = paleta) +
  theme(strip.text.x = element_text(face = "bold"),
        axis.title.y = element_text(face = "bold")) +
  theme(strip.text.x = element_text(face = "italic")) +
  theme(legend.text = element_text(face = "italic")) +
  theme(strip.background = element_rect(colour = "black", fill = "white")) +
  theme(text = element_text(size = 10))
## Warning in geom_bar(stat = "identity", pisition = "stack", res = 300): Ignoring
## unknown parameters: `pisition` and `res`
print(Rel_Abun_Phyl_Season)

#ggsave("../figures/Rel_Abun_Phyl_Season.jpeg", width=15, height=4.5, dpi=300)

PLOT RELATIVE ABUNDANCE AT THE GENUS LEVEL

Load data

metadata <- read.csv(file = "../data/metadata.csv", check.names = F)
asv_table <- read.csv("../data/feature_table.csv", check.names = F)
taxonomy <- read.csv("../data/taxonomy_03.csv", check.names = F) %>% unite(
  taxa, Kingdom:Species, remove = F, sep = ";")
asvtable_metadata <- asv_table %>% 
  inner_join(taxonomy)
## Joining with `by = join_by(OTUID)`

Relative abundance at phylum level (Sceloporus species)

Genus_01 <- asvtable_metadata %>% 
  group_by(Genus) %>% # But also for other taxonomic levels
  summarise_if(is.numeric, sum)
Genus_01 <- Genus_01 %>% 
  column_to_rownames(var = "Genus")
Genus.ra <- t(t(Genus_01)/colSums(Genus_01)*100)

rowMeans(Genus.ra) %>% 
  as.data.frame() %>% arrange(desc(.))
apply(Genus.ra,1,sd)
##                              Alistipes                            Bacteroides 
##                               2.748686                               8.221501 
##             Clostridia_vadinBB60_group                 Hafnia-Obesumbacterium 
##                               2.075043                               8.502137 
##                             Hungatella                      Lachnoclostridium 
##                               3.843127                               1.646522 
##                            Odoribacter                                 Others 
##                               3.661410                              12.190726 
##                        Parabacteroides                            Pseudomonas 
##                               3.747141                               5.199502 
##                              Roseburia  [Eubacterium]_coprostanoligenes_group 
##                               3.872776                               1.770494
Genus <- Genus.ra %>% t() %>% 
  as.data.frame()%>% 
  rownames_to_column(var = "SampleID") %>% 
  inner_join(metadata)
## Joining with `by = join_by(SampleID)`
Mean_Genus <- Genus %>% 
  group_by(Species_Scel) %>% 
  summarise_if(is.numeric, mean)
SD_Genus <- Genus %>% 
  group_by(Species_Scel) %>% 
  summarise_if(is.numeric, sd)
aggregate(Genus[ ,2:13], list(Genus$Species_Scel), mean)
aggregate(Genus[ ,2:13], list(Genus$Species_Scel), sd)
#write.table(Prom_Genus, file="./Genus_Pro_ST.txt", sep = "\t")

Filter genus taxonomic level

metadata <- read.csv(file = "../data/metadata.csv", header = TRUE, row.names = 1) 
asv_table <- read.csv(file = "../data/feature_table.csv", check.names = F) 
taxonomy <-  read.csv("../data/taxonomy_03.csv", check.names = F) %>% mutate_at(
  c("Genus"), str_replace,"g__", "")
lista <- rowMeans(Genus.ra) %>% as.data.frame() %>% 
  arrange(desc(.)) %>% 
  slice_head(n=14) %>% 
  rownames_to_column(var = "Genus") %>% 
  filter(!Genus =="g__") %>% 
  mutate_at(c("Genus"), str_replace,"g__", "")
list <- lista$Genus
lista01 <- read.csv(file = "../data/lista.csv", check.names = F)
list02 <- lista01$Genus
#write.table(lista, file="../data/lista.txt", sep = "\t")
taxonomy_filter <- taxonomy %>% 
  filter(Genus %in% list02)
taxonomy %>% inner_join(lista)
## Joining with `by = join_by(Genus)`
taxonomy_1 <- taxonomy_filter %>% 
  inner_join(asv_table, by =c(
  "OTUID"="OTUID")) %>% dplyr::select(1:8)

Load data again

asv_table_1 <- read.csv(file = "../data/feature_table.csv", header = TRUE,
                        row.names = 1) %>% 
  rownames_to_column(var = "OTUID") %>% 
  inner_join(taxonomy_1, by = "OTUID") %>% 
  dplyr::select(-97:-103) %>% 
  column_to_rownames(var = "OTUID")
taxo <- taxonomy_1 %>% column_to_rownames(var = "OTUID")

Create phyloseq objects

SAM <- sample_data(metadata)
TAX <- tax_table(as.matrix(taxo)) 
OTU <- otu_table(asv_table_1, taxa_are_rows=TRUE)  
physeq <- merge_phyloseq(OTU, TAX, SAM)

Convert to relative abundance

relative  = transform_sample_counts(physeq = physeq, function(OTU) OTU / sum(OTU))

Color palette

paleta <- c(brewer.pal(12, "Paired")[1:12], brewer.pal(8, "Dark2")[1:8])
print(paleta)
##  [1] "#A6CEE3" "#1F78B4" "#B2DF8A" "#33A02C" "#FB9A99" "#E31A1C" "#FDBF6F"
##  [8] "#FF7F00" "#CAB2D6" "#6A3D9A" "#FFFF99" "#B15928" "#1B9E77" "#D95F02"
## [15] "#7570B3" "#E7298A" "#66A61E" "#E6AB02" "#A6761D" "#666666"

Plot by Species and Season

Rel_Abun_Genus_Season <- plot_bar(physeq = relative, "Sample", fill = "Genus") + 
  facet_grid(~Species_Scel + Season, space = "free", scales = "free") +
  labs(y="Relative abundance") +
  geom_bar(stat = "identity", pisition = "stack", res = 300) +
  scale_fill_manual(values = paleta) +
  theme(strip.text.x = element_text(face = "bold"),
        axis.title.y = element_text(face = "bold")) +
  theme(strip.text.x = element_text(face = "italic")) +
  theme(legend.text = element_text(face = "italic")) +
  theme(strip.background = element_rect(colour = "black", fill = "white")) +
  theme(text = element_text(size = 10))
## Warning in geom_bar(stat = "identity", pisition = "stack", res = 300): Ignoring
## unknown parameters: `pisition` and `res`
print(Rel_Abun_Genus_Season)

#ggsave("../figures/Rel_Abun_Genus_Season.jpeg", width=15.0, height=6.5, dpi=400)

Core microbiota (50%)

aeneus_50 <- read.delim("../data/core_50_aeneus.tsv", sep = "\t", check.names = F) %>%
  rownames_to_column(var = "ids")

bicanthalis_50 <- read.delim("../data/core_50_bicanthalis.tsv", sep = "\t", check.names = F) %>%
  rownames_to_column(var = "ids")

grammicus_50 <- read.delim("../data/core_50_grammicus.tsv", sep = "\t", check.names = F) %>%
  rownames_to_column(var = "ids")

spinosus_50 <- read.delim("../data/core_50_spinosus.tsv", sep = "\t", check.names = F) %>%
  rownames_to_column(var = "ids") 

Plotting the Venn Diagramm

venn.plot_50 <- venn.diagram(
  x = list(S_grammicus = grammicus_50$OTUID,
           S_spinosus = spinosus_50$OTUID,
           S_aeneus = aeneus_50$OTUID,
           S_bicanthalis = bicanthalis_50$OTUID),
  category.names = c(
    expression(bold("S. grammicus")),
    expression(bold("S. spinosus")),
    expression(bold("S. aeneus")),
    expression(bold("S. bicanthalis"))),
  filename = "../figures/viendo_50.jpg",
  output = TRUE,
  height = 3000,
  width = 6000,
  resolution = 300,
  compression = "lzw",
  units = "px",
  lwd = 6,
  lty = "blank",
  fill = c("#56B4E9", "#009E73", "#999999", "#E69F00"),
  cex = 1.5,
  #fontface = "bold",
  fontface = "italic",
  fontfamily = "sans",
  cat.cex = 2,
  cat.fontface = "bold",
  cat.default.pos = "outer",
  cat.pos = c(-10, 25, -0.020, -0.050),
  cat.dist = c(0.3, 0.35, 0.099, 0.099),
  cat.fontfamily = "sans")

II. ALPHA DIVERSITY ANALYSES

Loading libraries

library(hilldiv)
library(ggplot2)
library(ggpubr)
library(reshape2)
library(tidyverse)
library(lme4)
library(car)
library(arm)
library(predictmeans)
library(AER)
library(visreg)
library(FSA)
library(rstatix)
library(sjPlot)
library(effects)
library(cowplot)

Calculating Hill numbers with hilldiv R package

otu_table <- read.csv("feature_table.csv", header = TRUE, row.names = 1)
q0 <- hill_div(otu_table, qvalue = 0)
q1 <- hill_div(otu_table, qvalue = 1)
q2 <- hill_div(otu_table, qvalue = 2)
q012 <- cbind(q0, q1, q2)
write.table(q012, file="../data/q012_hilldiv.txt", sep = "\t")

Preparing the data

Load files

richness_q012 <- read.csv("../data/Hill_numbers_q012.csv", header = TRUE) %>% 
  dplyr::select(SampleID, q0, q1, q2)
metadata <- read.csv("../data/metadata.csv", header = TRUE, check.names = F)
Micro_div <- richness_q012 %>% 
  inner_join(metadata, by = c("SampleID"="SampleID"))

Declare factors

SPECIES <- as.factor(Micro_div$Species_Scel) # Four Sceloporus lizard species
SEASON <- as.factor(Micro_div$Season) # Dry and Rainy
SEX <- as.factor(Micro_div$Sex) # Male and Female

Standardize continous independent variables

SVL <- rescale(Micro_div$SVL, binary.inputs = "center") # Snout–vent length measured in mm.
ELEVATION <- rescale(Micro_div$Elevation, binary.inputs = "center") # Taken as m a.s.l. 
SEQDEPTH <- rescale(Micro_div$SeqDepth, binary.inputs = "center")

Evaluate whether sex variable influence on gut microbiota diversity (q=1)

# Paired test (Wilcoxon test, q1) 

# Sceloporus aeneus
aeneus_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" & 
                            Sex == "Female", q1, drop = TRUE)
aeneus_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" & 
                             Sex == "Male", q1, drop = TRUE)
aeneus_FemMale_GutMic <- wilcox.test(x= aeneus_fem_GutMic, y= aeneus_male_GutMic)
aeneus_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  aeneus_fem_GutMic and aeneus_male_GutMic
## W = 49, p-value = 0.4779
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus bicanthalis
bica_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" & 
                          Sex == "Female", q1, drop = TRUE)
bica_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" & 
                           Sex == "Male", q1, drop = TRUE)
bica_FemMale_GutMic <- wilcox.test(x= bica_fem_GutMic, y= bica_male_GutMic)
bica_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  bica_fem_GutMic and bica_male_GutMic
## W = 69, p-value = 0.9095
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus grammicus
gram_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" & 
                          Sex == "Female", q1, drop = TRUE)
gram_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" & 
                           Sex == "Male", q1, drop = TRUE)
gram_FemMale_GutMic <- wilcox.test(x= gram_fem_GutMic, y= gram_male_GutMic)
gram_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  gram_fem_GutMic and gram_male_GutMic
## W = 109, p-value = 0.5676
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus spinosus
spi_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" & 
                         Sex == "Female", q1, drop = TRUE)
spi_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" & 
                          Sex == "Male", q1, drop = TRUE)
spi_FemMale_GutMic <- wilcox.test(x= spi_fem_GutMic, y= spi_male_GutMic)
spi_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  spi_fem_GutMic and spi_male_GutMic
## W = 55, p-value = 0.8596
## alternative hypothesis: true location shift is not equal to 0

Differences in snout vent length (SVL) among lizard species

set.seed(1234)
Data_SVL <- subset(Micro_div, select = c(Species_Scel, SVL))
kruskal.t <- Data_SVL %>% kruskal_test(SVL ~ Species_Scel)
kruskal.t
Post_hoc  <- Data_SVL %>% dunn_test(SVL ~ Species_Scel, p.adjust.method = "bonferroni") 
Post_hoc

Visualization: box plots with p-values

Post_hoc <- Post_hoc %>% 
  add_xy_position(x = "Species_Scel")
ggboxplot(Data_SVL, x = "Species_Scel", y = "SVL", fill = "Species_Scel") +
  xlab(element_blank())+
  scale_fill_manual(values = c("#56B4E9","#009E73", "#999999","#E69F00"))+
      theme_classic() +
   theme(legend.position = "right",
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        legend.title = element_blank()) +
  theme(legend.text = element_text(face = "italic"))+
  stat_pvalue_manual(Post_hoc, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(kruskal.t, detailed = TRUE),
    caption = get_pwc_label(Post_hoc))

Generalized Linear Model with quasi-Poisson distribution (q = 1)

M1 <- glm(q1 ~ SEASON*SPECIES+ELEVATION+SVL+SEQDEPTH,
           family = quasipoisson(link = "log"),
           data = Micro_div)
summary(M1)
## 
## Call:
## glm(formula = q1 ~ SEASON * SPECIES + ELEVATION + SVL + SEQDEPTH, 
##     family = quasipoisson(link = "log"), data = Micro_div)
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                4.262828   0.228362  18.667  < 2e-16
## SEASONRainy                                0.405260   0.151684   2.672 0.009061
## SPECIESSceloporus bicanthalis              0.989264   0.459589   2.152 0.034227
## SPECIESSceloporus grammicus                0.028572   0.168483   0.170 0.865747
## SPECIESSceloporus spinosus                 0.044405   0.217178   0.204 0.838486
## ELEVATION                                 -0.387007   0.399665  -0.968 0.335661
## SVL                                        0.002237   0.003986   0.561 0.576205
## SEQDEPTH                                   0.235518   0.065659   3.587 0.000561
## SEASONRainy:SPECIESSceloporus bicanthalis -0.641160   0.190178  -3.371 0.001132
## SEASONRainy:SPECIESSceloporus grammicus   -0.249934   0.196457  -1.272 0.206811
## SEASONRainy:SPECIESSceloporus spinosus    -0.239895   0.202098  -1.187 0.238567
##                                              
## (Intercept)                               ***
## SEASONRainy                               ** 
## SPECIESSceloporus bicanthalis             *  
## SPECIESSceloporus grammicus                  
## SPECIESSceloporus spinosus                   
## ELEVATION                                    
## SVL                                          
## SEQDEPTH                                  ***
## SEASONRainy:SPECIESSceloporus bicanthalis ** 
## SEASONRainy:SPECIESSceloporus grammicus      
## SEASONRainy:SPECIESSceloporus spinosus       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.08586)
## 
##     Null deviance: 1537.22  on 94  degrees of freedom
## Residual deviance:  999.94  on 84  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(M1)

M2 <- update(M1,. ~ . -SPECIES:SEASON)
anova(M1, M2, test = "F")
summary(M2)
## 
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + ELEVATION + SVL + SEQDEPTH, 
##     family = quasipoisson(link = "log"), data = Micro_div)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.462617   0.218130  20.458  < 2e-16 ***
## SEASONRainy                    0.114355   0.078245   1.462  0.14748    
## SPECIESSceloporus bicanthalis  0.652287   0.476593   1.369  0.17463    
## SPECIESSceloporus grammicus   -0.144594   0.116591  -1.240  0.21824    
## SPECIESSceloporus spinosus    -0.114180   0.186455  -0.612  0.54189    
## ELEVATION                     -0.467503   0.427006  -1.095  0.27661    
## SVL                            0.001759   0.004152   0.424  0.67284    
## SEQDEPTH                       0.205901   0.068543   3.004  0.00348 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.24543)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1136.0  on 87  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
M3 <- update(M2,. ~ . -SVL)
anova(M2, M3, test = "F")
summary(M3)
## 
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + ELEVATION + SEQDEPTH, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.53327    0.14023  32.327  < 2e-16 ***
## SEASONRainy                    0.11649    0.07772   1.499  0.13747    
## SPECIESSceloporus bicanthalis  0.66616    0.47539   1.401  0.16464    
## SPECIESSceloporus grammicus   -0.11718    0.09653  -1.214  0.22803    
## SPECIESSceloporus spinosus    -0.04897    0.10416  -0.470  0.63943    
## ELEVATION                     -0.48877    0.42406  -1.153  0.25220    
## SEQDEPTH                       0.20354    0.06796   2.995  0.00356 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.13221)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1138.2  on 88  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
M4 <- update(M3,. ~ . -ELEVATION) 
anova(M3, M4, test = "F")
summary(M4)
## 
## Call:
## glm(formula = q1 ~ SEASON + SPECIES + SEQDEPTH, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.66137    0.08481  54.964   <2e-16 ***
## SEASONRainy                    0.11993    0.07821   1.533   0.1287    
## SPECIESSceloporus bicanthalis  0.12860    0.09341   1.377   0.1721    
## SPECIESSceloporus grammicus   -0.11600    0.09699  -1.196   0.2349    
## SPECIESSceloporus spinosus    -0.01756    0.10109  -0.174   0.8625    
## SEQDEPTH                       0.20435    0.06833   2.991   0.0036 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.2577)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1157.3  on 89  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
confint(M4)
## Waiting for profiling to be done...
##                                     2.5 %     97.5 %
## (Intercept)                    4.49221597 4.82475124
## SEASONRainy                   -0.03336187 0.27330751
## SPECIESSceloporus bicanthalis -0.05390201 0.31251946
## SPECIESSceloporus grammicus   -0.30596309 0.07452011
## SPECIESSceloporus spinosus    -0.21612307 0.18046386
## SEQDEPTH                       0.06781350 0.33576942
M5 <- update(M4,. ~ . -SEQDEPTH) 
anova(M4, M5, test="F")
summary(M5)
## 
## Call:
## glm(formula = q1 ~ SEASON + SPECIES, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.60321    0.08660  53.155  < 2e-16 ***
## SEASONRainy                    0.22149    0.07245   3.057  0.00294 ** 
## SPECIESSceloporus bicanthalis  0.15745    0.09734   1.618  0.10927    
## SPECIESSceloporus grammicus   -0.11983    0.10129  -1.183  0.23993    
## SPECIESSceloporus spinosus    -0.01096    0.10574  -0.104  0.91766    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 13.41224)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1260.8  on 90  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Best Model
M6 <- glm(q1 ~ SEASON*SPECIES,
          family = quasipoisson(link = "log"),
          data = Micro_div)
anova(M6, test="F")
summary(M6)
## 
## Call:
## glm(formula = q1 ~ SEASON * SPECIES, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 4.3930     0.1316  33.371  < 2e-16
## SEASONRainy                                 0.5318     0.1561   3.406  0.00100
## SPECIESSceloporus bicanthalis               0.5414     0.1656   3.270  0.00154
## SPECIESSceloporus grammicus                 0.1005     0.1635   0.615  0.54032
## SPECIESSceloporus spinosus                  0.1593     0.1715   0.929  0.35550
## SEASONRainy:SPECIESSceloporus bicanthalis  -0.5796     0.2020  -2.869  0.00517
## SEASONRainy:SPECIESSceloporus grammicus    -0.3298     0.2065  -1.597  0.11392
## SEASONRainy:SPECIESSceloporus spinosus     -0.2364     0.2153  -1.098  0.27537
##                                              
## (Intercept)                               ***
## SEASONRainy                               ** 
## SPECIESSceloporus bicanthalis             ** 
## SPECIESSceloporus grammicus                  
## SPECIESSceloporus spinosus                   
## SEASONRainy:SPECIESSceloporus bicanthalis ** 
## SEASONRainy:SPECIESSceloporus grammicus      
## SEASONRainy:SPECIESSceloporus spinosus       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.61408)
## 
##     Null deviance: 1537.2  on 94  degrees of freedom
## Residual deviance: 1151.1  on 87  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
confint(M6)
## Waiting for profiling to be done...
##                                                2.5 %      97.5 %
## (Intercept)                                4.1233396  4.64031820
## SEASONRainy                                0.2317199  0.84503460
## SPECIESSceloporus bicanthalis              0.2207910  0.87128269
## SPECIESSceloporus grammicus               -0.2155594  0.42677045
## SPECIESSceloporus spinosus                -0.1742816  0.49967406
## SEASONRainy:SPECIESSceloporus bicanthalis -0.9789785 -0.18613994
## SEASONRainy:SPECIESSceloporus grammicus   -0.7385988  0.07179439
## SEASONRainy:SPECIESSceloporus spinosus    -0.6616867  0.18335432
coef(M6)
##                               (Intercept) 
##                                 4.3929537 
##                               SEASONRainy 
##                                 0.5317728 
##             SPECIESSceloporus bicanthalis 
##                                 0.5413804 
##               SPECIESSceloporus grammicus 
##                                 0.1005116 
##                SPECIESSceloporus spinosus 
##                                 0.1593282 
## SEASONRainy:SPECIESSceloporus bicanthalis 
##                                -0.5795744 
##   SEASONRainy:SPECIESSceloporus grammicus 
##                                -0.3297727 
##    SEASONRainy:SPECIESSceloporus spinosus 
##                                -0.2363617

Residual plot of the best model

Plot model assumption

plot_model(M6, type = "eff", terms = "SPECIES")

plot_model(M6, show.values = TRUE, value.offset = .3, width = 0.1,
           vline.color = "#E69F00")

residuals_q1 <- plot_model(M6, type = "eff", terms = c("SEASON","SPECIES"),
                           colors = c("#999999", "#E69F00", "#56B4E9", 
                                      "#009E73", "#F0E442")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  theme(axis.text.x = element_text(size = 11))
residuals_q1 

#ggsave("../figures/residuals_q1.jpeg", width = 5.0, height = 3.0, dpi = 300)

Contrasts (dry vs rainy), CI95%, p values

Tukey Test

library(emmeans)
emmeans(M6, specs = ~ SPECIES*SEASON) %>%
  contrast() %>%
  as.data.frame()
emm <- emmeans(M6, spec = ~ SPECIES*SEASON, 
               type = "response")
print(emm)
##  SPECIES                SEASON  rate    SE  df asymp.LCL asymp.UCL
##  Sceloporus aeneus      Dry     80.9 10.65 Inf      62.5       105
##  Sceloporus bicanthalis Dry    139.0 13.96 Inf     114.1       169
##  Sceloporus grammicus   Dry     89.4  8.67 Inf      74.0       108
##  Sceloporus spinosus    Dry     94.8 10.43 Inf      76.5       118
##  Sceloporus aeneus      Rainy  137.7 11.56 Inf     116.8       162
##  Sceloporus bicanthalis Rainy  132.5 10.56 Inf     113.3       155
##  Sceloporus grammicus   Rainy  109.4 10.31 Inf      91.0       132
##  Sceloporus spinosus    Rainy  127.4 12.68 Inf     104.9       155
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale
X <- contrast(emm, method = "pairwise")
confint(X)
M1.emm.s <- emmeans(M6, specs = ~ SPECIES*SEASON)
pairs(M1.emm.s, adjust = "tukey", infer=c(TRUE,TRUE))
##  contrast                                                  estimate    SE  df
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Dry        -0.54138 0.166 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Dry          -0.10051 0.164 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Dry           -0.15933 0.172 Inf
##  Sceloporus aeneus Dry - Sceloporus aeneus Rainy           -0.53177 0.156 Inf
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy      -0.49358 0.154 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Rainy        -0.30251 0.162 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Rainy         -0.45474 0.165 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Dry      0.44087 0.140 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Dry       0.38205 0.149 Inf
##  Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy       0.00961 0.131 Inf
##  Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy  0.04780 0.128 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy    0.23887 0.138 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy     0.08664 0.141 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Dry        -0.05882 0.147 Inf
##  Sceloporus grammicus Dry - Sceloporus aeneus Rainy        -0.43126 0.128 Inf
##  Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy   -0.39307 0.125 Inf
##  Sceloporus grammicus Dry - Sceloporus grammicus Rainy     -0.20200 0.135 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Rainy      -0.35423 0.139 Inf
##  Sceloporus spinosus Dry - Sceloporus aeneus Rainy         -0.37244 0.138 Inf
##  Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy    -0.33425 0.136 Inf
##  Sceloporus spinosus Dry - Sceloporus grammicus Rainy      -0.14318 0.145 Inf
##  Sceloporus spinosus Dry - Sceloporus spinosus Rainy       -0.29541 0.148 Inf
##  Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy     0.03819 0.116 Inf
##  Sceloporus aeneus Rainy - Sceloporus grammicus Rainy       0.22926 0.126 Inf
##  Sceloporus aeneus Rainy - Sceloporus spinosus Rainy        0.07703 0.130 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy  0.19107 0.123 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy   0.03884 0.127 Inf
##  Sceloporus grammicus Rainy - Sceloporus spinosus Rainy    -0.15223 0.137 Inf
##  asymp.LCL asymp.UCL z.ratio p.value
##    -1.0432   -0.0396  -3.270  0.0239
##    -0.5961    0.3950  -0.615  0.9987
##    -0.6792    0.3605  -0.929  0.9833
##    -1.0050   -0.0585  -3.406  0.0152
##    -0.9599   -0.0272  -3.208  0.0292
##    -0.7931    0.1880  -1.869  0.5723
##    -0.9548    0.0454  -2.756  0.1064
##     0.0178    0.8640   3.158  0.0341
##    -0.0693    0.8334   2.566  0.1686
##    -0.3871    0.4063   0.073  1.0000
##    -0.3407    0.4363   0.373  1.0000
##    -0.1784    0.6561   1.735  0.6640
##    -0.3418    0.5151   0.613  0.9987
##    -0.5032    0.3855  -0.401  0.9999
##    -0.8200   -0.0425  -3.362  0.0176
##    -0.7734   -0.0127  -3.132  0.0369
##    -0.6117    0.2077  -1.495  0.8107
##    -0.7753    0.0668  -2.550  0.1748
##    -0.7917    0.0469  -2.692  0.1248
##    -0.7458    0.0773  -2.462  0.2121
##    -0.5819    0.2956  -0.989  0.9761
##    -0.7448    0.1540  -1.992  0.4871
##    -0.3126    0.3890   0.330  1.0000
##    -0.1531    0.6116   1.817  0.6082
##    -0.3175    0.4716   0.592  0.9990
##    -0.1828    0.5649   1.549  0.7805
##    -0.3475    0.4251   0.305  1.0000
##    -0.5674    0.2629  -1.111  0.9546
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 8 estimates 
## P value adjustment: tukey method for comparing a family of 8 estimates
pwpm(M1.emm.s)
##                              Sceloporus aeneus Dry Sceloporus bicanthalis Dry
## Sceloporus aeneus Dry                       [4.39]                     0.0239
## Sceloporus bicanthalis Dry                -0.54138                     [4.93]
## Sceloporus grammicus Dry                  -0.10051                    0.44087
## Sceloporus spinosus Dry                   -0.15933                    0.38205
## Sceloporus aeneus Rainy                   -0.53177                    0.00961
## Sceloporus bicanthalis Rainy              -0.49358                    0.04780
## Sceloporus grammicus Rainy                -0.30251                    0.23887
## Sceloporus spinosus Rainy                 -0.45474                    0.08664
##                              Sceloporus grammicus Dry Sceloporus spinosus Dry
## Sceloporus aeneus Dry                          0.9987                  0.9833
## Sceloporus bicanthalis Dry                     0.0341                  0.1686
## Sceloporus grammicus Dry                       [4.49]                  0.9999
## Sceloporus spinosus Dry                      -0.05882                  [4.55]
## Sceloporus aeneus Rainy                      -0.43126                -0.37244
## Sceloporus bicanthalis Rainy                 -0.39307                -0.33425
## Sceloporus grammicus Rainy                   -0.20200                -0.14318
## Sceloporus spinosus Rainy                    -0.35423                -0.29541
##                              Sceloporus aeneus Rainy
## Sceloporus aeneus Dry                         0.0152
## Sceloporus bicanthalis Dry                    1.0000
## Sceloporus grammicus Dry                      0.0176
## Sceloporus spinosus Dry                       0.1248
## Sceloporus aeneus Rainy                       [4.92]
## Sceloporus bicanthalis Rainy                 0.03819
## Sceloporus grammicus Rainy                   0.22926
## Sceloporus spinosus Rainy                    0.07703
##                              Sceloporus bicanthalis Rainy
## Sceloporus aeneus Dry                              0.0292
## Sceloporus bicanthalis Dry                         1.0000
## Sceloporus grammicus Dry                           0.0369
## Sceloporus spinosus Dry                            0.2121
## Sceloporus aeneus Rainy                            1.0000
## Sceloporus bicanthalis Rainy                       [4.89]
## Sceloporus grammicus Rainy                        0.19107
## Sceloporus spinosus Rainy                         0.03884
##                              Sceloporus grammicus Rainy
## Sceloporus aeneus Dry                            0.5723
## Sceloporus bicanthalis Dry                       0.6640
## Sceloporus grammicus Dry                         0.8107
## Sceloporus spinosus Dry                          0.9761
## Sceloporus aeneus Rainy                          0.6082
## Sceloporus bicanthalis Rainy                     0.7805
## Sceloporus grammicus Rainy                       [4.70]
## Sceloporus spinosus Rainy                      -0.15223
##                              Sceloporus spinosus Rainy
## Sceloporus aeneus Dry                           0.1064
## Sceloporus bicanthalis Dry                      0.9987
## Sceloporus grammicus Dry                        0.1748
## Sceloporus spinosus Dry                         0.4871
## Sceloporus aeneus Rainy                         0.9990
## Sceloporus bicanthalis Rainy                    1.0000
## Sceloporus grammicus Rainy                      0.9546
## Sceloporus spinosus Rainy                       [4.85]
## 
## Row and column labels: SPECIES:SEASON
## Upper triangle: P values   adjust = "tukey"
## Diagonal: [Estimates] (emmean) 
## Lower triangle: Comparisons (estimate)   earlier vs. later
summary(M1.emm.s, type="response")
plot(M1.emm.s, comparisons=TRUE, type="response")

Evaluate whether sex variable influence on gut microbiota diversity (q=2)

# Paired test (Wilcoxon test, q2) 

# Sceloporus aeneus
aeneus_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" & 
                            Sex == "Female", q2, drop = TRUE)
aeneus_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus aeneus" & 
                             Sex == "Male", q2, drop = TRUE)
aeneus_FemMale_GutMic <- wilcox.test(x= aeneus_fem_GutMic, y= aeneus_male_GutMic)
aeneus_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  aeneus_fem_GutMic and aeneus_male_GutMic
## W = 39, p-value = 0.1713
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus bicanthalis
bica_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" & 
                          Sex == "Female", q2, drop = TRUE)
bica_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus bicanthalis" & 
                           Sex == "Male", q2, drop = TRUE)
bica_FemMale_GutMic <- wilcox.test(x= bica_fem_GutMic, y= bica_male_GutMic)
bica_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  bica_fem_GutMic and bica_male_GutMic
## W = 69, p-value = 0.9095
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus grammicus
gram_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" & 
                          Sex == "Female", q2, drop = TRUE)
gram_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus grammicus" & 
                           Sex == "Male", q2, drop = TRUE)
gram_FemMale_GutMic <- wilcox.test(x= gram_fem_GutMic, y= gram_male_GutMic)
gram_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  gram_fem_GutMic and gram_male_GutMic
## W = 107, p-value = 0.6313
## alternative hypothesis: true location shift is not equal to 0
# Sceloporus spinosus
spi_fem_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" & 
                         Sex == "Female", q2, drop = TRUE)
spi_male_GutMic <- subset(Micro_div, Species_Scel == "Sceloporus spinosus" & 
                          Sex == "Male", q2, drop = TRUE)
spi_FemMale_GutMic <- wilcox.test(x= spi_fem_GutMic, y= spi_male_GutMic)
spi_FemMale_GutMic
## 
##  Wilcoxon rank sum exact test
## 
## data:  spi_fem_GutMic and spi_male_GutMic
## W = 40, p-value = 0.4137
## alternative hypothesis: true location shift is not equal to 0

Generalized Linear Model with quasi-Poisson distribution (q = 2)

MODEL1 <- glm(q2 ~ SEASON*SPECIES+ELEVATION+SVL+SEQDEPTH,
            family = quasipoisson(link = "log"),
            data = Micro_div)
summary(MODEL1)
## 
## Call:
## glm(formula = q2 ~ SEASON * SPECIES + ELEVATION + SVL + SEQDEPTH, 
##     family = quasipoisson(link = "log"), data = Micro_div)
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                3.489e+00  3.263e-01  10.694
## SEASONRainy                                3.432e-01  2.143e-01   1.601
## SPECIESSceloporus bicanthalis              1.267e+00  6.553e-01   1.933
## SPECIESSceloporus grammicus                6.931e-02  2.358e-01   0.294
## SPECIESSceloporus spinosus                 1.270e-01  3.062e-01   0.415
## ELEVATION                                 -4.765e-01  5.739e-01  -0.830
## SVL                                       -2.387e-05  5.754e-03  -0.004
## SEQDEPTH                                   1.204e-01  1.009e-01   1.193
## SEASONRainy:SPECIESSceloporus bicanthalis -8.202e-01  2.652e-01  -3.093
## SEASONRainy:SPECIESSceloporus grammicus   -1.219e-01  2.754e-01  -0.442
## SEASONRainy:SPECIESSceloporus spinosus    -3.317e-01  2.909e-01  -1.140
##                                           Pr(>|t|)    
## (Intercept)                                < 2e-16 ***
## SEASONRainy                                0.11305    
## SPECIESSceloporus bicanthalis              0.05657 .  
## SPECIESSceloporus grammicus                0.76956    
## SPECIESSceloporus spinosus                 0.67944    
## ELEVATION                                  0.40867    
## SVL                                        0.99670    
## SEQDEPTH                                   0.23636    
## SEASONRainy:SPECIESSceloporus bicanthalis  0.00269 ** 
## SEASONRainy:SPECIESSceloporus grammicus    0.65931    
## SEASONRainy:SPECIESSceloporus spinosus     0.25746    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 9.481509)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  818.79  on 84  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(MODEL1)

MODEL2 <- update(MODEL1,. ~ . -SPECIES:SEASON)
anova(MODEL1, MODEL2, test = "F")
summary(MODEL2)
## 
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION + SVL + SEQDEPTH, 
##     family = quasipoisson(link = "log"), data = Micro_div)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.708501   0.314372  11.797   <2e-16 ***
## SEASONRainy                    0.025511   0.109991   0.232    0.817    
## SPECIESSceloporus bicanthalis  0.933649   0.676327   1.380    0.171    
## SPECIESSceloporus grammicus   -0.025499   0.165697  -0.154    0.878    
## SPECIESSceloporus spinosus    -0.051587   0.267991  -0.192    0.848    
## ELEVATION                     -0.626981   0.606572  -1.034    0.304    
## SVL                           -0.001225   0.005989  -0.205    0.838    
## SEQDEPTH                       0.041928   0.106778   0.393    0.696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.47637)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  934.33  on 87  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
MODEL3 <- update(MODEL2,. ~ . -SVL)
anova(MODEL2, MODEL3, test="F")
summary(MODEL3)
## 
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION + SEQDEPTH, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.65897    0.19913  18.375   <2e-16 ***
## SEASONRainy                    0.02422    0.10921   0.222    0.825    
## SPECIESSceloporus bicanthalis  0.92494    0.66933   1.382    0.171    
## SPECIESSceloporus grammicus   -0.04438    0.13685  -0.324    0.746    
## SPECIESSceloporus spinosus    -0.09673    0.15179  -0.637    0.526    
## ELEVATION                     -0.61281    0.59747  -1.026    0.308    
## SEQDEPTH                       0.04319    0.10604   0.407    0.685    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.36163)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  934.77  on 88  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
MODEL4 <- update(MODEL3,. ~ . -SEQDEPTH)
anova(MODEL3, MODEL4, test="F")
summary(MODEL4)
## 
## Call:
## glm(formula = q2 ~ SEASON + SPECIES + ELEVATION, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.64561    0.19533  18.664   <2e-16 ***
## SEASONRainy                    0.04488    0.09589   0.468    0.641    
## SPECIESSceloporus bicanthalis  0.93231    0.66505   1.402    0.164    
## SPECIESSceloporus grammicus   -0.04395    0.13615  -0.323    0.748    
## SPECIESSceloporus spinosus    -0.09462    0.15095  -0.627    0.532    
## ELEVATION                     -0.61383    0.59391  -1.034    0.304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.26101)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  936.46  on 89  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
MODEL5 <- update(MODEL4,. ~ . -ELEVATION) # Best model
anova(MODEL4, MODEL5, test="F")
summary(MODEL5)
## 
## Call:
## glm(formula = q2 ~ SEASON + SPECIES, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.80662    0.11637  32.712   <2e-16 ***
## SEASONRainy                    0.04940    0.09642   0.512   0.6097    
## SPECIESSceloporus bicanthalis  0.25823    0.13175   1.960   0.0531 .  
## SPECIESSceloporus grammicus   -0.04288    0.13674  -0.314   0.7546    
## SPECIESSceloporus spinosus    -0.05584    0.14700  -0.380   0.7049    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.3626)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  950.07  on 90  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Best Model
MODEL6 <- glm(q2 ~ SEASON*SPECIES,
              family = quasipoisson(link = "log"),
              data = Micro_div)
summary(MODEL6)
## 
## Call:
## glm(formula = q2 ~ SEASON * SPECIES, family = quasipoisson(link = "log"), 
##     data = Micro_div)
## 
## Coefficients:
##                                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                3.57869    0.17122  20.902  < 2e-16
## SEASONRainy                                0.40334    0.20706   1.948 0.054645
## SPECIESSceloporus bicanthalis              0.74146    0.20804   3.564 0.000596
## SPECIESSceloporus grammicus                0.08107    0.21339   0.380 0.704917
## SPECIESSceloporus spinosus                 0.15897    0.22310   0.713 0.478026
## SEASONRainy:SPECIESSceloporus bicanthalis -0.79089    0.26304  -3.007 0.003451
## SEASONRainy:SPECIESSceloporus grammicus   -0.14701    0.27125  -0.542 0.589223
## SEASONRainy:SPECIESSceloporus spinosus    -0.32729    0.29016  -1.128 0.262430
##                                              
## (Intercept)                               ***
## SEASONRainy                               .  
## SPECIESSceloporus bicanthalis             ***
## SPECIESSceloporus grammicus                  
## SPECIESSceloporus spinosus                   
## SEASONRainy:SPECIESSceloporus bicanthalis ** 
## SEASONRainy:SPECIESSceloporus grammicus      
## SEASONRainy:SPECIESSceloporus spinosus       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 9.452614)
## 
##     Null deviance: 1039.03  on 94  degrees of freedom
## Residual deviance:  840.53  on 87  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
coef(MODEL6)
##                               (Intercept) 
##                                3.57869413 
##                               SEASONRainy 
##                                0.40334451 
##             SPECIESSceloporus bicanthalis 
##                                0.74145835 
##               SPECIESSceloporus grammicus 
##                                0.08107416 
##                SPECIESSceloporus spinosus 
##                                0.15897437 
## SEASONRainy:SPECIESSceloporus bicanthalis 
##                               -0.79088870 
##   SEASONRainy:SPECIESSceloporus grammicus 
##                               -0.14700948 
##    SEASONRainy:SPECIESSceloporus spinosus 
##                               -0.32729337
confint(MODEL6)
## Waiting for profiling to be done...
##                                                  2.5 %     97.5 %
## (Intercept)                                3.223230360  3.8964819
## SEASONRainy                                0.006183591  0.8207875
## SPECIESSceloporus bicanthalis              0.342095601  1.1605800
## SPECIESSceloporus grammicus               -0.330307974  0.5093207
## SPECIESSceloporus spinosus                -0.274369614  0.6037278
## SEASONRainy:SPECIESSceloporus bicanthalis -1.313319633 -0.2802605
## SEASONRainy:SPECIESSceloporus grammicus   -0.685102597  0.3801804
## SEASONRainy:SPECIESSceloporus spinosus    -0.902353379  0.2372010

Residual plot of the best model

Plot model assumption

plot_model(MODEL6, type = "eff", terms = "SPECIES")

plot_model(MODEL6, show.values = TRUE, value.offset = .3, width = 0.1,
           vline.color = "#E69F00")

residuals_q2 <- plot_model(MODEL6, type = "eff", terms = c("SEASON","SPECIES"),
                           colors = c("#999999", "#E69F00", "#56B4E9", 
                                      "#009E73", "#F0E442")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  theme(axis.text.x = element_text(size = 11))
residuals_q2 

#ggsave("../figures/residuals_q2.jpeg", width = 5.0, height = 3.0, dpi = 300)

Contrasts (dry vs rainy), CI95%, p values

## Tukey Test

emmeans(MODEL6, specs = ~ SPECIES*SEASON) %>%
  contrast() %>%
  as.data.frame()
emm_q2 <- emmeans(MODEL6, spec = ~ SPECIES*SEASON, 
                  type = "response")
print(emm_q2)
##  SPECIES                SEASON rate   SE  df asymp.LCL asymp.UCL
##  Sceloporus aeneus      Dry    35.8 6.13 Inf      25.6      50.1
##  Sceloporus bicanthalis Dry    75.2 8.89 Inf      59.7      94.8
##  Sceloporus grammicus   Dry    38.9 4.95 Inf      30.3      49.9
##  Sceloporus spinosus    Dry    42.0 6.01 Inf      31.7      55.6
##  Sceloporus aeneus      Rainy  53.6 6.24 Inf      42.7      67.4
##  Sceloporus bicanthalis Rainy  51.0 5.67 Inf      41.1      63.5
##  Sceloporus grammicus   Rainy  50.2 6.04 Inf      39.7      63.6
##  Sceloporus spinosus    Rainy  45.3 6.55 Inf      34.1      60.1
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale
X_q2 <- contrast(emm_q2, method = "pairwise")
confint(X_q2)
M1.emm.s_q2 <- emmeans(MODEL6, specs = ~ SPECIES*SEASON)
pairs(M1.emm.s_q2, adjust = "tukey", infer=c(TRUE,TRUE))
##  contrast                                                  estimate    SE  df
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Dry         -0.7415 0.208 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Dry           -0.0811 0.213 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Dry            -0.1590 0.223 Inf
##  Sceloporus aeneus Dry - Sceloporus aeneus Rainy            -0.4033 0.207 Inf
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy       -0.3539 0.204 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Rainy         -0.3374 0.209 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Rainy          -0.2350 0.224 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Dry       0.6604 0.174 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Dry        0.5825 0.186 Inf
##  Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy        0.3381 0.166 Inf
##  Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy   0.3875 0.162 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy     0.4040 0.169 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy      0.5064 0.187 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Dry         -0.0779 0.192 Inf
##  Sceloporus grammicus Dry - Sceloporus aeneus Rainy         -0.3223 0.173 Inf
##  Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy    -0.2728 0.169 Inf
##  Sceloporus grammicus Dry - Sceloporus grammicus Rainy      -0.2563 0.175 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Rainy       -0.1540 0.193 Inf
##  Sceloporus spinosus Dry - Sceloporus aeneus Rainy          -0.2444 0.184 Inf
##  Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy     -0.1949 0.181 Inf
##  Sceloporus spinosus Dry - Sceloporus grammicus Rainy       -0.1784 0.187 Inf
##  Sceloporus spinosus Dry - Sceloporus spinosus Rainy        -0.0761 0.203 Inf
##  Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy      0.0494 0.161 Inf
##  Sceloporus aeneus Rainy - Sceloporus grammicus Rainy        0.0659 0.167 Inf
##  Sceloporus aeneus Rainy - Sceloporus spinosus Rainy         0.1683 0.186 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy   0.0165 0.164 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy    0.1189 0.182 Inf
##  Sceloporus grammicus Rainy - Sceloporus spinosus Rainy      0.1024 0.188 Inf
##  asymp.LCL asymp.UCL z.ratio p.value
##    -1.3720    -0.111  -3.564  0.0087
##    -0.7278     0.566  -0.380  0.9999
##    -0.8352     0.517  -0.713  0.9967
##    -1.0309     0.224  -1.948  0.5176
##    -0.9726     0.265  -1.734  0.6649
##    -0.9717     0.297  -1.612  0.7431
##    -0.9139     0.444  -1.049  0.9668
##     0.1338     1.187   3.801  0.0036
##     0.0201     1.145   3.139  0.0361
##    -0.1647     0.841   2.038  0.4560
##    -0.1041     0.879   2.389  0.2465
##    -0.1072     0.915   2.396  0.2433
##    -0.0592     1.072   2.714  0.1183
##    -0.6584     0.503  -0.407  0.9999
##    -0.8453     0.201  -1.868  0.5734
##    -0.7851     0.239  -1.614  0.7418
##    -0.7874     0.275  -1.463  0.8272
##    -0.7376     0.430  -0.800  0.9932
##    -0.8034     0.315  -1.325  0.8897
##    -0.7439     0.354  -1.076  0.9618
##    -0.7450     0.388  -0.955  0.9805
##    -0.6921     0.540  -0.374  1.0000
##    -0.4384     0.537   0.307  1.0000
##    -0.4416     0.573   0.394  0.9999
##    -0.3940     0.731   0.907  0.9855
##    -0.4799     0.513   0.101  1.0000
##    -0.4334     0.671   0.652  0.9981
##    -0.4674     0.672   0.545  0.9994
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 8 estimates 
## P value adjustment: tukey method for comparing a family of 8 estimates
pairs(M1.emm.s_q2, adjust = "tukey")
##  contrast                                                  estimate    SE  df
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Dry         -0.7415 0.208 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Dry           -0.0811 0.213 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Dry            -0.1590 0.223 Inf
##  Sceloporus aeneus Dry - Sceloporus aeneus Rainy            -0.4033 0.207 Inf
##  Sceloporus aeneus Dry - Sceloporus bicanthalis Rainy       -0.3539 0.204 Inf
##  Sceloporus aeneus Dry - Sceloporus grammicus Rainy         -0.3374 0.209 Inf
##  Sceloporus aeneus Dry - Sceloporus spinosus Rainy          -0.2350 0.224 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Dry       0.6604 0.174 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Dry        0.5825 0.186 Inf
##  Sceloporus bicanthalis Dry - Sceloporus aeneus Rainy        0.3381 0.166 Inf
##  Sceloporus bicanthalis Dry - Sceloporus bicanthalis Rainy   0.3875 0.162 Inf
##  Sceloporus bicanthalis Dry - Sceloporus grammicus Rainy     0.4040 0.169 Inf
##  Sceloporus bicanthalis Dry - Sceloporus spinosus Rainy      0.5064 0.187 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Dry         -0.0779 0.192 Inf
##  Sceloporus grammicus Dry - Sceloporus aeneus Rainy         -0.3223 0.173 Inf
##  Sceloporus grammicus Dry - Sceloporus bicanthalis Rainy    -0.2728 0.169 Inf
##  Sceloporus grammicus Dry - Sceloporus grammicus Rainy      -0.2563 0.175 Inf
##  Sceloporus grammicus Dry - Sceloporus spinosus Rainy       -0.1540 0.193 Inf
##  Sceloporus spinosus Dry - Sceloporus aeneus Rainy          -0.2444 0.184 Inf
##  Sceloporus spinosus Dry - Sceloporus bicanthalis Rainy     -0.1949 0.181 Inf
##  Sceloporus spinosus Dry - Sceloporus grammicus Rainy       -0.1784 0.187 Inf
##  Sceloporus spinosus Dry - Sceloporus spinosus Rainy        -0.0761 0.203 Inf
##  Sceloporus aeneus Rainy - Sceloporus bicanthalis Rainy      0.0494 0.161 Inf
##  Sceloporus aeneus Rainy - Sceloporus grammicus Rainy        0.0659 0.167 Inf
##  Sceloporus aeneus Rainy - Sceloporus spinosus Rainy         0.1683 0.186 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus grammicus Rainy   0.0165 0.164 Inf
##  Sceloporus bicanthalis Rainy - Sceloporus spinosus Rainy    0.1189 0.182 Inf
##  Sceloporus grammicus Rainy - Sceloporus spinosus Rainy      0.1024 0.188 Inf
##  z.ratio p.value
##   -3.564  0.0087
##   -0.380  0.9999
##   -0.713  0.9967
##   -1.948  0.5176
##   -1.734  0.6649
##   -1.612  0.7431
##   -1.049  0.9668
##    3.801  0.0036
##    3.139  0.0361
##    2.038  0.4560
##    2.389  0.2465
##    2.396  0.2433
##    2.714  0.1183
##   -0.407  0.9999
##   -1.868  0.5734
##   -1.614  0.7418
##   -1.463  0.8272
##   -0.800  0.9932
##   -1.325  0.8897
##   -1.076  0.9618
##   -0.955  0.9805
##   -0.374  1.0000
##    0.307  1.0000
##    0.394  0.9999
##    0.907  0.9855
##    0.101  1.0000
##    0.652  0.9981
##    0.545  0.9994
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 8 estimates
pwpm(M1.emm.s_q2)
##                              Sceloporus aeneus Dry Sceloporus bicanthalis Dry
## Sceloporus aeneus Dry                       [3.58]                     0.0087
## Sceloporus bicanthalis Dry                 -0.7415                     [4.32]
## Sceloporus grammicus Dry                   -0.0811                     0.6604
## Sceloporus spinosus Dry                    -0.1590                     0.5825
## Sceloporus aeneus Rainy                    -0.4033                     0.3381
## Sceloporus bicanthalis Rainy               -0.3539                     0.3875
## Sceloporus grammicus Rainy                 -0.3374                     0.4040
## Sceloporus spinosus Rainy                  -0.2350                     0.5064
##                              Sceloporus grammicus Dry Sceloporus spinosus Dry
## Sceloporus aeneus Dry                          0.9999                  0.9967
## Sceloporus bicanthalis Dry                     0.0036                  0.0361
## Sceloporus grammicus Dry                       [3.66]                  0.9999
## Sceloporus spinosus Dry                       -0.0779                  [3.74]
## Sceloporus aeneus Rainy                       -0.3223                 -0.2444
## Sceloporus bicanthalis Rainy                  -0.2728                 -0.1949
## Sceloporus grammicus Rainy                    -0.2563                 -0.1784
## Sceloporus spinosus Rainy                     -0.1540                 -0.0761
##                              Sceloporus aeneus Rainy
## Sceloporus aeneus Dry                         0.5176
## Sceloporus bicanthalis Dry                    0.4560
## Sceloporus grammicus Dry                      0.5734
## Sceloporus spinosus Dry                       0.8897
## Sceloporus aeneus Rainy                       [3.98]
## Sceloporus bicanthalis Rainy                  0.0494
## Sceloporus grammicus Rainy                    0.0659
## Sceloporus spinosus Rainy                     0.1683
##                              Sceloporus bicanthalis Rainy
## Sceloporus aeneus Dry                              0.6649
## Sceloporus bicanthalis Dry                         0.2465
## Sceloporus grammicus Dry                           0.7418
## Sceloporus spinosus Dry                            0.9618
## Sceloporus aeneus Rainy                            1.0000
## Sceloporus bicanthalis Rainy                       [3.93]
## Sceloporus grammicus Rainy                         0.0165
## Sceloporus spinosus Rainy                          0.1189
##                              Sceloporus grammicus Rainy
## Sceloporus aeneus Dry                            0.7431
## Sceloporus bicanthalis Dry                       0.2433
## Sceloporus grammicus Dry                         0.8272
## Sceloporus spinosus Dry                          0.9805
## Sceloporus aeneus Rainy                          0.9999
## Sceloporus bicanthalis Rainy                     1.0000
## Sceloporus grammicus Rainy                       [3.92]
## Sceloporus spinosus Rainy                        0.1024
##                              Sceloporus spinosus Rainy
## Sceloporus aeneus Dry                           0.9668
## Sceloporus bicanthalis Dry                      0.1183
## Sceloporus grammicus Dry                        0.9932
## Sceloporus spinosus Dry                         1.0000
## Sceloporus aeneus Rainy                         0.9855
## Sceloporus bicanthalis Rainy                    0.9981
## Sceloporus grammicus Rainy                      0.9994
## Sceloporus spinosus Rainy                       [3.81]
## 
## Row and column labels: SPECIES:SEASON
## Upper triangle: P values   adjust = "tukey"
## Diagonal: [Estimates] (emmean) 
## Lower triangle: Comparisons (estimate)   earlier vs. later
summary(M1.emm.s_q2, type="response")
plot(M1.emm.s_q2, comparisons=TRUE, type="response")

Relation between alpha taxonomic microbiota and diet

Micro_div$Ind <- substr(Micro_div$SampleID, 8, nchar(Micro_div$SampleID) - 0)

diet_phylo_div <- read.csv("../data/diet_Hill_numbers_q012.csv", header = TRUE) 
# Generar una nueva columna con la información de otra columna, excluyendo los últimos 4 caracteres
diet_phylo_div$Ind <- substr(diet_phylo_div$SampleID, 8, nchar(diet_phylo_div$SampleID) - 4)
#names(diet_phylo_div)[names(diet_phylo_div) %in% c("q0", "q1", "q2")] <- c("d.p.0", "d.p.1", #"d.p.2")
divs <- Micro_div %>% 
  inner_join(diet_phylo_div, by = c("Ind"="Ind"))

#write.table(divs, file="../data/full_divs.txt", sep = "\t")

##Plot

taxq1_tq0 <- ggscatter(divs, x = "tq0", y = "q1", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), xlab= "Taxonomic richness of diet (number of taxonomic groups)", ylab="Gut bacterial diversity (effective number of common ASVs)",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
   )+
  theme(legend.title = element_blank(), legend.position = "none")


taxq2_tq0 <- ggscatter(divs, x = "tq0", y = "q2", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Taxonomic richness of diet (number of taxonomic groups)", ylab="Gut bacterial diversity (effective number of dominant ASVs)",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
   )+
  theme(legend.title = element_blank(), legend.position = "none")


phyq1_pq0 <- ggscatter(divs, x = "pq0", y = "q1", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73") , 
   xlab= " Phylogenetic richness of diet (number of phylogenetic lineages)", ylab="Gut bacterial diversity (effective number of common ASVs)",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
   )+
  theme(legend.title = element_blank(), legend.position = "none")


phyq2_pq0 <- ggscatter(divs, x = "pq0", y = "q2", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Phylogenetic richness of diet (number of phylogenetic lineages)", 
   ylab="Gut bacterial diversity (effective number of dominant ASVs)",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
   )+
  theme(legend.title = element_blank(), legend.position = "none")


taxq1_tq0.sp <- ggscatter(divs, x = "tq0", y = "q1", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Taxonomic richness of diet (number of taxonomic groups)", 
   ylab="Gut bacterial diversity (effective number of common ASVs)",
   facet.by = "Species_Scel",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", 
                         label.x = 3, label.sep = "\n")
   )+
  theme(legend.position = "none",
        strip.text = element_text(face = "italic"))



taxq2_tq0.sp <- ggscatter(divs, x = "tq0", y = "q2", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Taxonomic richness of diet (number of taxonomic groups)", ylab="Gut bacterial diversity (effective number of dominant ASVs)",
   facet.by = "Species_Scel" ,
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
   )+
  theme(legend.position = "none",
        strip.text = element_text(face = "italic"))


phyq1_pq0.sp <- ggscatter(divs, x = "pq0", y = "q1", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73") , 
   xlab= " Phylogenetic richness of diet (number of phylogenetic lineages)", ylab="Gut bacterial diversity (effective number of common ASVs)",
   facet.by = "Species_Scel",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
   )+
  theme(legend.position = "none",
        strip.text = element_text(face = "italic"))


phyq2_pq0.sp <- ggscatter(divs, x = "pq0", y = "q2", color = "Species_Scel",
   palette = c("#999999", "#E69F00", "#56B4E9", "#009E73"), 
   xlab= "Phylogenetic richness of diet (number of phylogenetic lineages)", ylab="Gut bacterial diversity (effective number of dominant ASVs)", 
   facet.by = "Species_Scel" ,
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE, # Add confidence interval
   cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
   cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
   )+
  theme(legend.position = "none",
        strip.text = element_text(face = "italic"))

taxq1_tq0

taxq1_tq0.sp

taxq2_tq0

taxq2_tq0.sp

phyq1_pq0

phyq1_pq0.sp

phyq2_pq0

phyq2_pq0.sp

diet.vs.micro1<- plot_grid(taxq1_tq0, taxq1_tq0.sp, 
                          labels = "AUTO", ncol = 1)
#ggsave("../figures/dietvsmicro1.jpeg", width=8, height=11, dpi=300)
diet.vs.micro2<- plot_grid(taxq2_tq0, taxq2_tq0.sp, 
                          labels = c("E","F"), ncol = 1)
#ggsave("../figures/dietvsmicro2.jpeg", width=8, height=11, dpi=300)
diet.vs.micro3<- plot_grid(phyq1_pq0, phyq1_pq0.sp, 
                          labels = c("C", "D"), ncol = 1)
#ggsave("../figures/dietvsmicro3.jpeg", width=8, height=11, dpi=300)
diet.vs.micro4<- plot_grid(phyq2_pq0, phyq2_pq0.sp, 
                          labels = c("G", "H"), ncol = 1)
#ggsave("../figures/dietvsmicro4.jpeg", width=8, height=11, dpi=300)

Effect of diet on gut microbiota

loading files

richness_microb <- read.csv("../data/Hill_numbers_q012.csv", header = TRUE) %>% 
  dplyr::select(SampleID, q0, q1, q2)
metadata <- read.csv("../data/Metadata_model_DietMic.csv", header = TRUE, check.names = F)
Microb_data <- richness_microb %>% inner_join(metadata, by = c("SampleID"="SampleID"))

Check distribution of the dependent variable (diversity q=0, Hill numbers)

shapiro.test(Microb_data$q0) # Data are not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  Microb_data$q0
## W = 0.92649, p-value = 0.0009384
shapiro.test(Microb_data$q1) # Normal data
## 
##  Shapiro-Wilk normality test
## 
## data:  Microb_data$q1
## W = 0.97555, p-value = 0.2334
hist(Microb_data$q0, col = "darkgreen") # q0

hist(log(Microb_data$q0), col = "orange") # log transformation

hist(Microb_data$q1, col = "gray") # q1

hist(Microb_data$q2, col = "blue") # q2

bartlett.test(q0 ~ Species_Scel, data = Microb_data) 
## 
##  Bartlett test of homogeneity of variances
## 
## data:  q0 by Species_Scel
## Bartlett's K-squared = 9.0931, df = 3, p-value = 0.02808
leveneTest(q0 ~ Species_Scel, data = Microb_data)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

Plot variances among lizard species (q0)

boxplot(q0 ~ Species_Scel, data = Microb_data,
        main = "Differences in variance among lizard species",
        xlab = "Species",
        ylab = " Hill numbers (q0)",
        col = "#4682B433",
        border = "black")

###Preparing the data

#Declare as factors
SPECIES <- as.factor(Microb_data$Species_Scel) # Four Sceloporus lizard species
SEASON <- as.factor(Microb_data$Season) # Dry and Rainy
SEX <- as.factor(Microb_data$Sex) # Male and Female

#Standardize continous independent variables
SVL <- rescale(Microb_data$SVL, binary.inputs = "center") # Snout–vent length measured in mm.
ELEVATION <- rescale(Microb_data$Elevation, binary.inputs = "center") # Taken as m a.s.l. 
SEQDEPTH <- rescale(Microb_data$SeqDepth, binary.inputs = "center")

Linear model with Poisson distribution (taxonomic q = 0)

With Microbial Diversity (q1)

## UPDATE FUNCTION (RUN ANOVA)
# 68 samples (diet / microbiota)
m1 <- glm(q1 ~ Diet*SPECIES+SEASON,
          family = quasipoisson(link = "log"),
          data = Microb_data)
summary(m1)
## 
## Call:
## glm(formula = q1 ~ Diet * SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         4.6562151  0.1865713  24.957  < 2e-16 ***
## Diet                               -0.0004683  0.0196703  -0.024  0.98109    
## SPECIESSceloporus bicanthalis       0.0595147  0.2216353   0.269  0.78930    
## SPECIESSceloporus grammicus        -0.1948473  0.2453897  -0.794  0.43059    
## SPECIESSceloporus spinosus          0.1751032  0.2349154   0.745  0.45921    
## SEASONRainy                         0.2689010  0.0857259   3.137  0.00274 ** 
## Diet:SPECIESSceloporus bicanthalis  0.0018106  0.0258806   0.070  0.94448    
## Diet:SPECIESSceloporus grammicus   -0.0084125  0.0312962  -0.269  0.78909    
## Diet:SPECIESSceloporus spinosus    -0.0373851  0.0318303  -1.175  0.24525    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.57198)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 593.29  on 55  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(m1)

m2 <- update(m1,. ~ . -Diet:SPECIES)
anova(m1, m2, test="F")
summary(m2)
## 
## Call:
## glm(formula = q1 ~ Diet + SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.704855   0.130358  36.092  < 2e-16 ***
## Diet                          -0.007877   0.010173  -0.774  0.44187    
## SPECIESSceloporus bicanthalis  0.061506   0.106646   0.577  0.56636    
## SPECIESSceloporus grammicus   -0.256866   0.111237  -2.309  0.02452 *  
## SPECIESSceloporus spinosus    -0.047488   0.116892  -0.406  0.68605    
## SEASONRainy                    0.283586   0.080238   3.534  0.00081 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.42645)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 614.21  on 58  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m3 <- update(m2,. ~ . -Diet)
anova(m2, m3, test="F")
summary(m3) # Best Model
## 
## Call:
## glm(formula = q1 ~ SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.63657    0.09649  48.053  < 2e-16 ***
## SPECIESSceloporus bicanthalis  0.07079    0.10567   0.670 0.505539    
## SPECIESSceloporus grammicus   -0.24165    0.10918  -2.213 0.030761 *  
## SPECIESSceloporus spinosus    -0.02901    0.11425  -0.254 0.800452    
## SEASONRainy                    0.29232    0.07918   3.692 0.000488 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.37193)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 620.50  on 59  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m4 <- update(m3,. ~ . -SPECIES)
anova(m3, m4, test="F")
summary(m4)
## 
## Call:
## glm(formula = q1 ~ SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.55181    0.06272  72.574  < 2e-16 ***
## SEASONRainy  0.33820    0.08005   4.225 7.99e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.18809)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 728.20  on 62  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m5 <- update(m3,. ~ . -SEASON)
anova(m3, m5, test="F")
summary(m5)
## 
## Call:
## glm(formula = q1 ~ SPECIES, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.82632    0.08770  55.033   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.09084    0.11575   0.785   0.4356    
## SPECIESSceloporus grammicus   -0.29554    0.11870  -2.490   0.0156 *  
## SPECIESSceloporus spinosus    -0.08290    0.12430  -0.667   0.5074    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 12.47322)
## 
##     Null deviance: 932.15  on 63  degrees of freedom
## Residual deviance: 764.28  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Plot model assumption
plot_model(m3, type = "eff", terms = "SPECIES",
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme(legend.text = element_text(face = "italic"))

plot_model(m3, type = "eff", terms = c("SPECIES", "SEASON"),
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  theme(axis.text.x = element_text(size = 11))

With Microbial Diversity (q2)

model1 <- glm(q2 ~ Diet*SPECIES+SEASON,
              family = quasipoisson(link = "log"),
              data = Microb_data)
summary(model1)
## 
## Call:
## glm(formula = q2 ~ Diet * SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         3.62441    0.25085  14.449   <2e-16 ***
## Diet                                0.01764    0.02568   0.687   0.4950    
## SPECIESSceloporus bicanthalis       0.34411    0.29428   1.169   0.2473    
## SPECIESSceloporus grammicus         0.08006    0.32076   0.250   0.8038    
## SPECIESSceloporus spinosus          0.45688    0.31449   1.453   0.1520    
## SEASONRainy                         0.11588    0.10904   1.063   0.2925    
## Diet:SPECIESSceloporus bicanthalis -0.01057    0.03311  -0.319   0.7507    
## Diet:SPECIESSceloporus grammicus   -0.03106    0.04019  -0.773   0.4429    
## Diet:SPECIESSceloporus spinosus    -0.07581    0.04204  -1.803   0.0768 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.255196)
## 
##     Null deviance: 528.4  on 63  degrees of freedom
## Residual deviance: 401.7  on 55  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(model1)

model2 <- update(model1,. ~ . -Diet:SPECIES)
anova(model1, model2, test="F")
summary(model2)
## 
## Call:
## glm(formula = q2 ~ Diet + SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.79639    0.17233  22.030   <2e-16 ***
## Diet                          -0.00507    0.01324  -0.383   0.7030    
## SPECIESSceloporus bicanthalis  0.24216    0.14164   1.710   0.0927 .  
## SPECIESSceloporus grammicus   -0.15074    0.14887  -1.013   0.3155    
## SPECIESSceloporus spinosus    -0.02013    0.15835  -0.127   0.8993    
## SEASONRainy                    0.13451    0.10404   1.293   0.2012    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.364622)
## 
##     Null deviance: 528.40  on 63  degrees of freedom
## Residual deviance: 430.43  on 58  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model3 <- update(model2,. ~ . -Diet)
anova(model2, model3, test="F")
summary(model3)
## 
## Call:
## glm(formula = q2 ~ SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.752299   0.127845  29.350   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.247615   0.139789   1.771   0.0817 .  
## SPECIESSceloporus grammicus   -0.140309   0.145194  -0.966   0.3378    
## SPECIESSceloporus spinosus    -0.008757   0.154400  -0.057   0.9550    
## SEASONRainy                    0.140365   0.102034   1.376   0.1741    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.248433)
## 
##     Null deviance: 528.40  on 63  degrees of freedom
## Residual deviance: 431.52  on 59  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model4 <- update(model3,. ~ . -SEASON)
anova(model3, model4, test="F") # Best model
summary(model4)
## 
## Call:
## glm(formula = q2 ~ SPECIES, family = quasipoisson(link = "log"), 
##     data = Microb_data)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.84098    0.10960  35.046   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.25751    0.13984   1.841   0.0705 .  
## SPECIESSceloporus grammicus   -0.16641    0.14420  -1.154   0.2531    
## SPECIESSceloporus spinosus    -0.03486    0.15350  -0.227   0.8211    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.272215)
## 
##     Null deviance: 528.40  on 63  degrees of freedom
## Residual deviance: 445.32  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model5 <- update(model4,. ~ . -SPECIES)
anova(model4, model5, test="F")
summary(model5)
## 
## Call:
## glm(formula = q2 ~ 1, family = quasipoisson(link = "log"), data = Microb_data)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.8565     0.0515   74.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 8.029905)
## 
##     Null deviance: 528.4  on 63  degrees of freedom
## Residual deviance: 528.4  on 63  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Plot model assumption
plot_model(model3, type = "eff", terms = "SPECIES",
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme(legend.text = element_text(face = "italic"))

Linear models with Poisson distribution (Phylogenetic q = 0)

With Microbial Diversity (q1)

divs$phyDiet<- divs$pq0
#Declare as factors
SPECIES <- as.factor(divs$Species_Scel) # Four Sceloporus lizard species
SEASON <- as.factor(divs$Season) # Dry and Rainy
SEX <- as.factor(divs$Sex) # Male and Female

#Standardize continous independent variables
SVL <- rescale(divs$SVL, binary.inputs = "center") # Snout–vent length measured in mm.
ELEVATION <- rescale(divs$Elevation, binary.inputs = "center") # Taken as m a.s.l. 
SEQDEPTH <- rescale(divs$SeqDepth, binary.inputs = "center")

## UPDATE FUNCTION (RUN ANOVA)
# 68 samples (diet / microbiota)
m1 <- glm(q1 ~ phyDiet*SPECIES+SEASON,
          family = quasipoisson(link = "log"),
          data = divs)
summary(m1)
## 
## Call:
## glm(formula = q1 ~ phyDiet * SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            4.661408   0.206567  22.566  < 2e-16 ***
## phyDiet                               -0.005446   0.069609  -0.078  0.93791    
## SPECIESSceloporus bicanthalis         -0.042498   0.236501  -0.180  0.85803    
## SPECIESSceloporus grammicus           -0.262152   0.243544  -1.076  0.28628    
## SPECIESSceloporus spinosus            -0.013347   0.247408  -0.054  0.95717    
## SEASONRainy                            0.276381   0.088922   3.108  0.00294 ** 
## phyDiet:SPECIESSceloporus bicanthalis  0.050503   0.084691   0.596  0.55332    
## phyDiet:SPECIESSceloporus grammicus   -0.004840   0.086901  -0.056  0.95577    
## phyDiet:SPECIESSceloporus spinosus    -0.014879   0.082267  -0.181  0.85711    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.85698)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 619.92  on 57  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(m1)

m2 <- update(m1,. ~ . -phyDiet:SPECIES)
anova(m1, m2, test="F")
summary(m2)
## 
## Call:
## glm(formula = q1 ~ phyDiet + SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.651537   0.133638  34.807  < 2e-16 ***
## phyDiet                        0.001562   0.026733   0.058  0.95360    
## SPECIESSceloporus bicanthalis  0.073228   0.106615   0.687  0.49483    
## SPECIESSceloporus grammicus   -0.274588   0.113016  -2.430  0.01812 *  
## SPECIESSceloporus spinosus    -0.062852   0.109585  -0.574  0.56842    
## SEASONRainy                    0.264348   0.086185   3.067  0.00324 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.49499)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 633.01  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m3 <- update(m2,. ~ . -phyDiet)
anova(m2, m3, test="F")
summary(m3) # Best Model
## 
## Call:
## glm(formula = q1 ~ SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.65694    0.09566  48.684  < 2e-16 ***
## SPECIESSceloporus bicanthalis  0.07275    0.10543   0.690  0.49280    
## SPECIESSceloporus grammicus   -0.27534    0.11136  -2.472  0.01622 *  
## SPECIESSceloporus spinosus    -0.06223    0.10818  -0.575  0.56720    
## SEASONRainy                    0.26231    0.07814   3.357  0.00136 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 10.32373)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 633.05  on 61  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m4 <- update(m3,. ~ . -SPECIES)
anova(m3, m4, test="F")
summary(m4)
## 
## Call:
## glm(formula = q1 ~ SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.55181    0.06312  72.115  < 2e-16 ***
## SEASONRainy  0.32085    0.07995   4.013  0.00016 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.33081)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 758.28  on 64  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
m5 <- update(m3,. ~ . -SEASON)
anova(m3, m5, test="F")
summary(m5)
## 
## Call:
## glm(formula = q1 ~ SPECIES, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.82632    0.08573  56.296  < 2e-16 ***
## SPECIESSceloporus bicanthalis  0.09084    0.11315   0.803  0.42514    
## SPECIESSceloporus grammicus   -0.33142    0.11834  -2.801  0.00679 ** 
## SPECIESSceloporus spinosus    -0.08424    0.11604  -0.726  0.47057    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.91967)
## 
##     Null deviance: 944.88  on 65  degrees of freedom
## Residual deviance: 751.44  on 62  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Plot model assumption
plot_model(m3, type = "eff", terms = "SPECIES",
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme(legend.text = element_text(face = "italic"))

plot_model(m3, type = "eff", terms = c("SPECIES", "SEASON"),
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  theme(axis.text.x = element_text(size = 11))

With Microbial Diversity (q2)

model1 <- glm(q2 ~ phyDiet*SPECIES+SEASON,
              family = quasipoisson(link = "log"),
              data = divs)
summary(model1)
## 
## Call:
## glm(formula = q2 ~ phyDiet * SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            3.878108   0.280118  13.845   <2e-16 ***
## phyDiet                               -0.042399   0.096811  -0.438    0.663    
## SPECIESSceloporus bicanthalis          0.005418   0.319419   0.017    0.987    
## SPECIESSceloporus grammicus           -0.258721   0.327267  -0.791    0.432    
## SPECIESSceloporus spinosus             0.045784   0.335106   0.137    0.892    
## SEASONRainy                            0.112807   0.115877   0.974    0.334    
## phyDiet:SPECIESSceloporus bicanthalis  0.101063   0.114513   0.883    0.381    
## phyDiet:SPECIESSceloporus grammicus    0.034825   0.117544   0.296    0.768    
## phyDiet:SPECIESSceloporus spinosus    -0.026253   0.113567  -0.231    0.818    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.702367)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 447.64  on 57  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
plot(model1)

model2 <- update(model1,. ~ . -phyDiet:SPECIES)
anova(model1, model2, test="F")
summary(model2)
## 
## Call:
## glm(formula = q2 ~ phyDiet + SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.82070    0.17848  21.407   <2e-16 ***
## phyDiet                       -0.01355    0.03571  -0.379   0.7058    
## SPECIESSceloporus bicanthalis  0.24687    0.14407   1.714   0.0918 .  
## SPECIESSceloporus grammicus   -0.17655    0.15307  -1.153   0.2533    
## SPECIESSceloporus spinosus    -0.04872    0.15269  -0.319   0.7508    
## SEASONRainy                    0.08802    0.11442   0.769   0.4448    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.670608)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 467.11  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model3 <- update(model2,. ~ . -phyDiet)
anova(model2, model3, test="F")
summary(model3)
## 
## Call:
## glm(formula = q2 ~ SPECIES + SEASON, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.77434    0.12992  29.052   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.24998    0.14286   1.750   0.0852 .  
## SPECIESSceloporus grammicus   -0.17020    0.15117  -1.126   0.2646    
## SPECIESSceloporus spinosus    -0.05479    0.15090  -0.363   0.7178    
## SEASONRainy                    0.10615    0.10332   1.027   0.3083    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.57)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 468.22  on 61  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model4 <- update(model3,. ~ . -SEASON)
anova(model3, model4, test="F") # Best model
summary(model4)
## 
## Call:
## glm(formula = q2 ~ SPECIES, family = quasipoisson(link = "log"), 
##     data = divs)
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.84098    0.11157  34.427   <2e-16 ***
## SPECIESSceloporus bicanthalis  0.25751    0.14236   1.809   0.0753 .  
## SPECIESSceloporus grammicus   -0.19303    0.14920  -1.294   0.2006    
## SPECIESSceloporus spinosus    -0.06383    0.15031  -0.425   0.6725    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 7.536168)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 476.25  on 62  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
model5 <- update(model4,. ~ . -SPECIES)
anova(model4, model5, test="F")
summary(model5)
## 
## Call:
## glm(formula = q2 ~ 1, family = quasipoisson(link = "log"), data = divs)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.84306    0.05214    73.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 8.3741)
## 
##     Null deviance: 568.65  on 65  degrees of freedom
## Residual deviance: 568.65  on 65  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Plot model assumption
plot_model(model3, type = "eff", terms = "SPECIES",
           colors = c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442")) +
  theme(legend.text = element_text(face = "italic"))

III. BETA DIVERSITY ANALYSES

Install hilldiv2

install.packages("devtools")
library(devtools)
install_github("anttonalberdi/hilldiv2")

Loading libraries

library(hilldiv2)
library(tidyverse)
library(vegan)
library(phytools)
library(pairwiseAdonis)
library(reshape2)
library(ggpubr)
library(rstatix)
# Load files
otutable <- read.csv("../data/feature_table.csv", row.names = 1)
Metadata <- read.csv("../data/metadata_beta.csv")

Calculate beta diversity

hill_pair_dis <- hillpair(data = otutable, q=c(1,2))
## The estimated time for calculating the 4465 pairwise combinations is 5 seconds.
names(hill_pair_dis)
## [1] "q1S" "q1C" "q1U" "q1V" "q2S" "q2C" "q2U" "q2V"
head(hill_pair_dis$q1S)
## [1] 0.6500844 0.6843612 0.7557878 0.1858421 0.5226403 0.6284240
  • S (Jaccard-type turnover): it quantifies the normalised species turnover rate in a sample relative to the total pool of samples (i.e., gamma diversity).
  • V (Sørensen-type turnover): it quantifies the normalised species turnover rate in a sample relative to one sample (i.e., alpha diversity).
  • U (Jaccard-type overlap-complement): it quantifies the proportion of non-shared species in the total pool of samples. Therefore, this metric quantifies dissimilarity from the perspective of the pool of samples.
  • C (Sørensen-type overlap-complement): it quantifies the effective average proportion of non-shared OTUs/ASVs/MAGs in samples. Therefore, this metric quantifies dissimilarity from the perspective of a single sample.

Taxonomic beta diversity at order q1

hill_pair_dis_nmds <- hill_pair_dis$q1S %>% #select the distance object based on dissimilarity metric S
  metaMDS(.,trymax = 500, k = 2, verbosity = FALSE) %>%
  vegan:::scores.metaMDS() %>%
  as_tibble(., rownames = "sample")

hill_pair_dis_nmds <- hill_pair_dis_nmds %>%
  left_join(Metadata, by = join_by(sample == SampleID)) %>%
  group_by(Species, Season) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup()

Plot ordination q1

NMDS_Plot_q1 <- ggplot(hill_pair_dis_nmds, aes(x = NMDS1, y = NMDS2, 
                                               color = Species, shape = Season)) +
  geom_point(size = 2) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.2) +
  theme(legend.position = "right",
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank()) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  scale_colour_manual(values = c("#999999","#E69F00","#56B4E9","#009E73","#F0E442"))
print(NMDS_Plot_q1)

#ggsave("../figures/NMDS_Plot_q1.jpeg", width = 8.0, height = 4.0, dpi = 300)

Run PERMANONA

set.seed(123)
dist <- hill_pair_dis$q1S
perm_q1 <- adonis2(dist ~ Species*Season+Elevation+SVL+SeqDepth,
                   data = Metadata, permutations = 999)
print(perm_q1)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist ~ Species * Season + Elevation + SVL + SeqDepth, data = Metadata, permutations = 999)
##                Df SumOfSqs      R2       F Pr(>F)    
## Species         3   6.6636 0.22329 10.2239  0.001 ***
## Season          1   1.7738 0.05944  8.1645  0.001 ***
## Elevation       1   0.2071 0.00694  0.9531  0.566    
## SVL             1   0.3386 0.01134  1.5583  0.053 .  
## SeqDepth        1   0.2866 0.00960  1.3191  0.118    
## Species:Season  3   2.3235 0.07786  3.5648  0.001 ***
## Residual       84  18.2495 0.61153                   
## Total          94  29.8426 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Run pairwise perMANOVA

Metadata$species_season <- paste(Metadata$Species, Metadata$Season, sep = "_")
pw_adonis_q1<-pairwise.adonis(dist,factors=Metadata$species_season)

intrasp_compar <- pw_adonis_q1[grep("^([^_]+)_Dry vs \\1_Rainy", pw_adonis_q1$pairs), ]
print(intrasp_compar)
##                                                         pairs Df SumsOfSqs
## 5      Sceloporus grammicus_Dry vs Sceloporus grammicus_Rainy  1 1.2919702
## 10       Sceloporus spinosus_Dry vs Sceloporus spinosus_Rainy  1 0.9619991
## 17           Sceloporus aeneus_Dry vs Sceloporus aeneus_Rainy  1 0.9720931
## 22 Sceloporus bicanthalis_Dry vs Sceloporus bicanthalis_Rainy  1 0.9324159
##     F.Model        R2 p.value p.adjusted sig
## 5  5.919983 0.1854632   0.001      0.028   .
## 10 3.704441 0.1631593   0.001      0.028   .
## 17 4.424882 0.1811629   0.001      0.028   .
## 22 5.104701 0.1883327   0.001      0.028   .

Plot of intraspecific distances

compar_intrasp<- c("Sceloporus grammicus_Dry_vs_Sceloporus grammicus_Rainy", "Sceloporus aeneus_Dry_vs_Sceloporus aeneus_Rainy", "Sceloporus bicanthalis_Dry_vs_Sceloporus bicanthalis_Rainy", "Sceloporus spinosus_Dry_vs_Sceloporus spinosus_Rainy")

dist.q1 <- 1-(hill_pair_dis$q1S) %>% as.matrix()

#convert matrix into df
dist.q1[lower.tri(dist.q1)] <- NA
dist.q1.df <- as.matrix(dist.q1) %>% melt(
  varnames = c("site1", "site2"))%>% drop_na() %>% filter(!value==0)

#Join beta diversity vakues with metadata
dist.q1.meta<- dist.q1.df %>% 
  inner_join(Metadata, by = c("site1"="SampleID")) %>% 
  inner_join(Metadata, by = c("site2"="SampleID"))

#Subset to obtain specific comparisons 
dist.q1.meta.filt <- dist.q1.meta %>% unite(
  "compar_intraspp", c("species_season.x", "species_season.y"), sep="_vs_", remove=F) %>% filter(compar_intraspp %in% compar_intrasp) 

# Plotting

#Setting the order of the species
dist.q1.meta.filt$Species.order<- factor(dist.q1.meta.filt$Species.x,
                                         levels = c("Sceloporus aeneus", "Sceloporus bicanthalis", "Sceloporus grammicus", "Sceloporus spinosus"))

# Kruskal wallis test
set.seed(1234)
kruskal.t.q1 <- dist.q1.meta.filt %>% kruskal_test(value ~ Species.order)
Post_hoc.q1  <- dist.q1.meta.filt %>% dunn_test(value ~ Species.order, p.adjust.method = "bonferroni") 
Post_hoc.q1 <- Post_hoc.q1 %>% 
  add_xy_position(x = "Species.order")

# Boxplot
q1_intrasp<-dist.q1.meta.filt %>% 
  ggboxplot(x = "Species.order", y="value", fill = "Species.order")+
  ylab("Seasonal turnover of common ASVs (q=1)")+
  xlab(element_blank())+
  scale_fill_manual(values = c("#999999","#E69F00","#56B4E9","#009E73","#F0E442"))+
      theme_classic() +
   theme(legend.position = "right",
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        legend.title = element_blank()) +
  theme(legend.text = element_text(face = "italic"))+
  stat_pvalue_manual(Post_hoc.q1, 
                     hide.ns = TRUE, 
                     step.increase = 0.2) +
  labs(
    subtitle = get_test_label(kruskal.t.q1, 
                              detailed = TRUE),
    caption = get_pwc_label(Post_hoc.q1))+
  scale_y_continuous(limits = c(0,1))

  
q1_intrasp

#ggsave("../figures/beta_div_intrasp_q1.jpeg", width = 5.0, height = 4.0, dpi = 300)

Run BETADISPER analysis

disp.species = betadisper(dist, Metadata$Species)
anova(disp.species)
permutest(disp.species)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)   
## Groups     3 0.06544 0.0218120 4.1684    999  0.008 **
## Residuals 91 0.47618 0.0052328                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permutest(disp.species, pairwise = TRUE, permutations = 999)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
## Groups     3 0.06544 0.0218120 4.1684    999  0.013 *
## Residuals 91 0.47618 0.0052328                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##                        Sceloporus aeneus Sceloporus bicanthalis
## Sceloporus aeneus                                    0.17700000
## Sceloporus bicanthalis        0.16097306                       
## Sceloporus grammicus          0.66564466             0.01087924
## Sceloporus spinosus           0.11780891             0.00011976
##                        Sceloporus grammicus Sceloporus spinosus
## Sceloporus aeneus                0.66800000               0.122
## Sceloporus bicanthalis           0.01700000               0.001
## Sceloporus grammicus                                      0.091
## Sceloporus spinosus              0.09456538
boxplot(disp.species)

Taxonomic beta diversity at order q2

hill_pair_dis_nmds_q2 <- hill_pair_dis$q2S %>% #select the distance object based on dissimilarity metric S
  metaMDS(.,trymax = 500, k = 2, verbosity = FALSE) %>%
  vegan:::scores.metaMDS() %>%
  as_tibble(., rownames = "sample")

hill_pair_dis_nmds_2 <- hill_pair_dis_nmds_q2 %>%
  left_join(Metadata, by = join_by(sample == SampleID)) %>%
  group_by(Species, Season) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup()

Plot ordination q2

NMDS_Plot_q2 <- ggplot(hill_pair_dis_nmds_2, aes(x = NMDS1, y = NMDS2, 
                                                  color = Species, shape = Season)) +
  geom_point(size = 2) +
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.2) +
  theme(legend.position = "right",
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank()) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic")) +
  scale_colour_manual(values = c("#999999","#E69F00","#56B4E9","#009E73","#F0E442"))
print(NMDS_Plot_q2)

#ggsave("../figures/NMDS_Plot_q2.jpeg", width = 8.0, height = 4.0, dpi = 300)

Run PERMANOVA

set.seed(123)
dist2 <- hill_pair_dis$q2S
perm_q2 <- adonis2(dist2 ~ Species*Season+Elevation+SVL+SeqDepth,
                   data = Metadata, permutations = 999)
print(perm_q2)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist2 ~ Species * Season + Elevation + SVL + SeqDepth, data = Metadata, permutations = 999)
##                Df SumOfSqs      R2      F Pr(>F)    
## Species         3    6.665 0.19570 8.7137  0.001 ***
## Season          1    2.129 0.06251 8.3503  0.001 ***
## Elevation       1    0.226 0.00664 0.8867  0.659    
## SVL             1    0.393 0.01155 1.5431  0.067 .  
## SeqDepth        1    0.294 0.00865 1.1548  0.270    
## Species:Season  3    2.933 0.08612 3.8346  0.001 ***
## Residual       84   21.418 0.62884                  
## Total          94   34.060 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Run pairwise perMANOVA

pw_adonis_q2<-pairwise.adonis(dist2,factors=Metadata$species_season)
intrasp_compar_q2 <- pw_adonis_q2[grep("^([^_]+)_Dry vs \\1_Rainy", pw_adonis_q2$pairs), ]
print(intrasp_compar_q2)
##                                                         pairs Df SumsOfSqs
## 5      Sceloporus grammicus_Dry vs Sceloporus grammicus_Rainy  1 1.3213474
## 10       Sceloporus spinosus_Dry vs Sceloporus spinosus_Rainy  1 0.9123456
## 17           Sceloporus aeneus_Dry vs Sceloporus aeneus_Rainy  1 1.2367957
## 22 Sceloporus bicanthalis_Dry vs Sceloporus bicanthalis_Rainy  1 1.6336642
##     F.Model        R2 p.value p.adjusted sig
## 5  5.294296 0.1691777   0.001      0.028   .
## 10 2.860417 0.1308492   0.001      0.028   .
## 17 4.641605 0.1883646   0.001      0.028   .
## 22 8.146037 0.2702192   0.001      0.028   .

Plot of intraspecific distances

compar_intrasp<- c("Sceloporus grammicus_Dry_vs_Sceloporus grammicus_Rainy", "Sceloporus aeneus_Dry_vs_Sceloporus aeneus_Rainy", "Sceloporus bicanthalis_Dry_vs_Sceloporus bicanthalis_Rainy", "Sceloporus spinosus_Dry_vs_Sceloporus spinosus_Rainy")

dist.q2 <- 1-(hill_pair_dis$q2S) %>% as.matrix()

#convert matrix into df
dist.q2[lower.tri(dist.q2)] <- NA
dist.q2.df <- as.matrix(dist.q2) %>% melt(
  varnames = c("site1", "site2"))%>% drop_na() %>% filter(!value==0)

#Join beta diversity vakues with metadata
dist.q2.meta<- dist.q2.df %>% 
  inner_join(Metadata, by = c("site1"="SampleID")) %>% 
  inner_join(Metadata, by = c("site2"="SampleID"))

#Subset to obtain specific comparisons 
dist.q2.meta.filt <- dist.q2.meta %>% unite(
  "compar_intraspp", c("species_season.x", "species_season.y"), sep="_vs_", remove=F) %>% filter(compar_intraspp %in% compar_intrasp) 

# Plotting

#Setting the order of the species
dist.q2.meta.filt$Species.order<- factor(dist.q2.meta.filt$Species.x,
                                         levels = c("Sceloporus aeneus", "Sceloporus bicanthalis", "Sceloporus grammicus", "Sceloporus spinosus"))

# Kruskal wallis test
set.seed(1234)
kruskal.t.q2 <- dist.q2.meta.filt %>% kruskal_test(value ~ Species.order)
Post_hoc.q2  <- dist.q2.meta.filt %>% dunn_test(value ~ Species.order, p.adjust.method = "bonferroni") 
Post_hoc.q2 <- Post_hoc.q2 %>% 
  add_xy_position(x = "Species.order")

q2_intrasp<-dist.q2.meta.filt %>% 
 ggboxplot(x = "Species.order", y="value", fill = "Species.order")+
  ylab("Seasonal turnover of dominant ASVs (q=2)")+
  xlab(element_blank())+
  scale_fill_manual(values = c("#999999","#E69F00","#56B4E9","#009E73","#F0E442"))+
      theme_classic() +
   theme(legend.position = "right",
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        legend.title = element_blank()) +
  theme(legend.text = element_text(face = "italic"))+
  stat_pvalue_manual(Post_hoc.q2, 
                     hide.ns = TRUE,
                     step.increase = 0.1) +
  labs(
    subtitle = get_test_label(kruskal.t.q2, 
                              detailed = TRUE),
    caption = get_pwc_label(Post_hoc.q2))+
  scale_y_continuous(limits = c(0,1))
  
q2_intrasp

#ggsave("../figures/beta_div_intrasp_q2.jpeg", width = 5.0, height = 4.0, dpi = 300)

Run BETADISPER analysis

disp.species_q2 = betadisper(dist2, Metadata$Species)
anova(disp.species_q2)
permutest(disp.species_q2)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
## Groups     3 0.07338 0.0244593 2.5963    999  0.043 *
## Residuals 91 0.85729 0.0094208                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permutest(disp.species_q2, pairwise = TRUE, permutations = 999)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
## Groups     3 0.07338 0.0244593 2.5963    999  0.067 .
## Residuals 91 0.85729 0.0094208                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##                        Sceloporus aeneus Sceloporus bicanthalis
## Sceloporus aeneus                                     0.2130000
## Sceloporus bicanthalis         0.2201605                       
## Sceloporus grammicus           0.6839133              0.2468309
## Sceloporus spinosus            0.2294011              0.0059616
##                        Sceloporus grammicus Sceloporus spinosus
## Sceloporus aeneus                 0.6960000               0.232
## Sceloporus bicanthalis            0.2560000               0.007
## Sceloporus grammicus                                      0.030
## Sceloporus spinosus               0.0297543
boxplot(disp.species_q2)

IV. MANTEL Y PROCRUSTES ANALYSIS

Load libraries

library(vegan)
library(ggplot2)
library(cowplot)
library(corrr)
library(tidyverse)
library(ape)
library(hilldiv2)

Loading files

This analysis is based on 68 samples obtained from both datasets (diet-microbiota)

metadata <- read.csv("../data/Procrustes_metadata.csv", header = TRUE, check.names = FALSE, 
                     row.names = 1)
diet <- read.csv("../data/Procrustes_Diet.csv", header = TRUE, check.names = FALSE, 
                 row.names = 1) %>% t()
microbiota <- read.csv("../data/Procrustes_Microbiota.csv", header = TRUE, 
                       check.names = FALSE, row.names = 1)
#microbiota <- microbiota %>% filter(rowSums(across(where(is.numeric)))!=0) %>% t()

Calculate distances matrices

# Distance matrix of bacterial microbiota
jaccard_mic_1 <- hillpair(data = microbiota, q = 1) # distance matrix at order q1
microbiota_jaccard_1 <- as.dist(jaccard_mic_1$S)
jaccard_mic_2 <- hillpair(data = microbiota, q = 2) # distance matrix at order q1
microbiota_jaccard_2 <- as.dist(jaccard_mic_2$S)

# Distance matrix of diet
diet_jaccard <- vegdist(diet, method = "jaccard")

Mantel test

## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = diet_jaccard, ydis = microbiota_jaccard_1, method = "pearson",      permutations = 999) 
## 
## Mantel statistic r: 0.206 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0349 0.0473 0.0560 0.0687 
## Permutation: free
## Number of permutations: 999
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = diet_jaccard, ydis = microbiota_jaccard_2, method = "pearson",      permutations = 999) 
## 
## Mantel statistic r: 0.1948 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0387 0.0511 0.0610 0.0764 
## Permutation: free
## Number of permutations: 999

Make pcoas

#PcoA of diet
dietPCoA <- as.data.frame(cmdscale(diet_jaccard))
plot(dietPCoA)

# PCoA of microbiota q=1
microbiotaPCoA_1 <- as.data.frame(cmdscale(microbiota_jaccard_1))
plot(microbiotaPCoA_1)

# PCoA of microbiota q=2
microbiotaPCoA_2 <- as.data.frame(cmdscale(microbiota_jaccard_2))
plot(microbiotaPCoA_2)

Procrustes analysis of all species

Procrustes analysis of all species at order q=1

procrust <- procrustes(X = dietPCoA, Y = microbiotaPCoA_1, scale=TRUE,
                       symmetric = TRUE)
pro_test <- protest(dietPCoA, microbiotaPCoA_1, permutations = 9999)
pro_test
## 
## Call:
## protest(X = dietPCoA, Y = microbiotaPCoA_1, permutations = 9999) 
## 
## Procrustes Sum of Squares (m12 squared):        0.8268 
## Correlation in a symmetric Procrustes rotation: 0.4161 
## Significance:  1e-04 
## 
## Permutation: free
## Number of permutations: 9999
plot(pro_test)

eigen <- sqrt(procrust$svd$d)
percent_var <- signif(eigen/sum(eigen), 4)*100

beta_pro <- data.frame(procrust$X)
head(beta_pro)
trans_pro <- data.frame(procrust$Yrot)
head(trans_pro)
beta_pro$UserName <- rownames(beta_pro)
beta_pro$type <- "Diet"
beta_pro$seasons<-metadata[,5]
beta_pro$species<-metadata[,1]
head(beta_pro)
trans_pro$UsarName <- rownames(trans_pro)
trans_pro$type <- "Microbiota"
trans_pro$seasons=metadata[,5]
trans_pro$species=metadata[,1]
head(trans_pro)
colnames(trans_pro) <- colnames(beta_pro)
toplot <- data.frame(rbind(beta_pro,trans_pro))
pval <- signif(pro_test$signif, 1)

#col1 = rgb(250/255,60/255,60/255)
#col4 = rgb(0/255,200/255,200/255)
#col8 = rgb(160/255,0/255,200/255)
#col10 = rgb(0/255,160/255,255/255)
# "#999999","#E69F00","#56B4E9", "#009E73","#F0E442

Graph diet and microbiota

procrustes_new <- ggplot(toplot) +
  geom_point(size = 1.8, alpha = 0.75, aes(x = V1, y = V2, color = species, 
                                         shape = type)) + 
  scale_color_manual(values = c("#999999","#E69F00","#56B4E9",
                                "#009E73")) +
  theme_classic() +
  scale_x_continuous(limits = c(-0.1, 0.3)) +
  scale_y_continuous(limits = c(-0.2, 0.2)) +
  geom_line(aes(x= V1, y = V2, group = UserName), 
            col = "darkgrey", alpha = 0.6, 
            size = 0.2) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 10, colour = "black"),
        legend.position = "right",
        axis.text = element_text(size = 10, colour = "black"),
        axis.title = element_text(size = 12, colour = "black"),
        aspect.ratio = 1) +
  theme(legend.text = element_text(face = "italic")) +
  guides(color = guide_legend(ncol = 1)) +
  annotate("text", x = 0.2, y = -0.2, 
           label = paste0("p-value=", pval), size = 4) +
  xlab(paste0("PCoA 1 (",percent_var[1],"%)")) +
  ylab(paste0("PCoA 2 (",percent_var[2],"%)"))
print(procrustes_new)

#ggsave("../figures/procrustes_diet_q1.jpeg", width = 8.0, height = 4.0, dpi = 300)

Procrustes analysis of all species at order q=2

procrust <- procrustes(X = dietPCoA, Y = microbiotaPCoA_2, scale=TRUE,
                       symmetric = TRUE)
pro_test <- protest(dietPCoA, microbiotaPCoA_2, permutations = 9999)
pro_test
## 
## Call:
## protest(X = dietPCoA, Y = microbiotaPCoA_2, permutations = 9999) 
## 
## Procrustes Sum of Squares (m12 squared):        0.8177 
## Correlation in a symmetric Procrustes rotation: 0.4269 
## Significance:  1e-04 
## 
## Permutation: free
## Number of permutations: 9999
plot(pro_test)

eigen <- sqrt(procrust$svd$d)
percent_var <- signif(eigen/sum(eigen), 4)*100

beta_pro <- data.frame(procrust$X)
head(beta_pro)
trans_pro <- data.frame(procrust$Yrot)
head(trans_pro)
beta_pro$UserName <- rownames(beta_pro)
beta_pro$type <- "Diet"
beta_pro$seasons<-metadata[,5]
beta_pro$species<-metadata[,1]
head(beta_pro)
trans_pro$UsarName <- rownames(trans_pro)
trans_pro$type <- "Microbiota"
trans_pro$seasons=metadata[,5]
trans_pro$species=metadata[,1]
head(trans_pro)
colnames(trans_pro) <- colnames(beta_pro)
toplot <- data.frame(rbind(beta_pro,trans_pro))
pval <- signif(pro_test$signif, 1)

Graph diet and microbiota

procrustes_new <- ggplot(toplot) +
  geom_point(size = 1.8, alpha = 0.75, aes(x = V1, y = V2, color = species, 
                                         shape = type)) + 
  scale_color_manual(values = c("#999999","#E69F00","#56B4E9",
                                "#009E73")) +
  theme_classic() +
  scale_x_continuous(limits = c(-0.1, 0.3)) +
  scale_y_continuous(limits = c(-0.2, 0.2)) +
  geom_line(aes(x= V1, y = V2, group = UserName), 
            col = "darkgrey", alpha = 0.6, 
            size = 0.2) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 10, colour = "black"),
        legend.position = "right",
        axis.text = element_text(size = 10, colour = "black"),
        axis.title = element_text(size = 12, colour = "black"),
        aspect.ratio = 1) +
  theme(legend.text = element_text(face = "italic")) +
  guides(color = guide_legend(ncol = 1)) +
  annotate("text", x = 0.2, y = -0.2, 
           label = paste0("p-value=", pval), size = 4) +
  xlab(paste0("PCoA 1 (",percent_var[1],"%)")) +
  ylab(paste0("PCoA 2 (",percent_var[2],"%)"))
print(procrustes_new)

#ggsave("../figures/procrustes_diet_q2.jpeg", width = 8.0, height = 4.0, dpi = 300)

Procrustes analysis by lizard species

Procrsutes by lizard species were done at q=1 as similar patterns were observed in the general procrustes at q=2.

Procrustes of Sceloporus aeneus

#loading data
aeneus_metadata <- subset(metadata, Species_Scel == "Sceloporus aeneus")
diet_aeneus <- read.csv("../data/Proc_Diet_aeneus.csv", header = TRUE, check.names = FALSE, 
                        row.names = 1) %>% t()
microbiota_aeneus <- read.csv("../data/Proc_Mic_aeneus.csv", header = TRUE, 
                              check.names = FALSE, row.names = 1)

# Distance matrix
jaccard_mic_a <- hillpair(data = microbiota_aeneus, q = 1)
microbiota_jaccard_a <- as.dist(jaccard_mic_a$S)

# Absence/presence data frame - Jaccard distance
diet_jaccard_a <- vegdist(diet_aeneus, method = "jaccard")
diet_jaccard_a <- as.dist(diet_jaccard_a)

# Mantel test

set.seed(123)
mantel(diet_jaccard_a, microbiota_jaccard_a, method = "pearson",
       permutations = 999)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = diet_jaccard_a, ydis = microbiota_jaccard_a, method = "pearson",      permutations = 999) 
## 
## Mantel statistic r: 0.3412 
##       Significance: 0.031 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.196 0.283 0.348 0.399 
## Permutation: free
## Number of permutations: 999
# Make pcoas

dietPCoA_a <- as.data.frame(cmdscale(diet_jaccard_a))
plot(dietPCoA_a)

microbiotaPCoA_a <- as.data.frame(cmdscale(microbiota_jaccard_a))
plot(microbiotaPCoA_a)

# Procrustes analysis

procrust_a <- procrustes(X = dietPCoA_a, Y = microbiotaPCoA_a, scale=TRUE,
                         symmetric = TRUE)
pro_test_a <- protest(dietPCoA_a, microbiotaPCoA_a, permutations = 9999)
pro_test_a
## 
## Call:
## protest(X = dietPCoA_a, Y = microbiotaPCoA_a, permutations = 9999) 
## 
## Procrustes Sum of Squares (m12 squared):        0.7374 
## Correlation in a symmetric Procrustes rotation: 0.5125 
## Significance:  0.0734 
## 
## Permutation: free
## Number of permutations: 9999
plot(pro_test_a)

eigen_a <- sqrt(procrust_a$svd$d)
percent_var_a <- signif(eigen_a/sum(eigen_a), 4)*100

beta_pro_a <- data.frame(procrust_a$X)
head(beta_pro_a)
trans_pro_a <- data.frame(procrust_a$Yrot)
head(trans_pro_a)
beta_pro_a$UserName <- rownames(beta_pro_a)
beta_pro_a$type <- "Diet"
beta_pro_a$seasons <- aeneus_metadata[,5]
beta_pro_a$species <- aeneus_metadata[,1]
head(beta_pro_a)
trans_pro_a$UsarName <- rownames(trans_pro_a)
trans_pro_a$type <- "Microbiota"
trans_pro_a$seasons <- aeneus_metadata[,5]
trans_pro_a$species <- aeneus_metadata[,1]
head(trans_pro_a)
colnames(trans_pro_a) <- colnames(beta_pro_a)
toplot_a <- data.frame(rbind(beta_pro_a, trans_pro_a))
pval_a <- signif(pro_test_a$signif, 1)


# Graph diet and microbiota

procreustes_aeneus_new <- ggplot(toplot_a) +
  geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species, 
                                         shape = type)) + 
  scale_color_manual(values = c("#999999")) +
  theme_classic() +
  scale_x_continuous(limits = c(-0.2, 0.4)) +
  scale_y_continuous(limits = c(-0.2, 0.4)) +
  geom_line(aes(x = V1, y = V2, group = UserName), 
            col = "darkgrey", alpha = 0.6, 
            size = 0.2) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 10, colour = "black"),
        legend.position = "right",
        axis.text = element_text(size = 10, colour = "black"),
        axis.title = element_text(size = 12, colour = "black"),
        aspect.ratio = 1) +
  theme(legend.text = element_text(face = "italic")) +
  guides(color = guide_legend(ncol = 1)) +
  #annotate("text", x = 0.07, y = -0.13, label = paste0("p-value=", pval_a), 
  #         size = 4) +
  xlab(paste0("PCoA 1 (",percent_var_a[1],"%)")) +
  ylab(paste0("PCoA 2 (",percent_var_a[2],"%)"))
print(procreustes_aeneus_new)

#ggsave("../figures/procreustes_aeneus_new.jpeg", width = 8.0, height = 4.0, dpi = 300)

Procrustes of Sceloporus bicanthalis

#Loading data
bica_metadata <- subset(metadata, Species_Scel == "Sceloporus bicanthalis")
diet_bica <- read.csv("../data/Proc_Diet_bica.csv", header = TRUE, check.names = FALSE, 
                      row.names = 1) %>% t()
microbiota_bica <- read.csv("../data/Proc_Mic_bica.csv", header = TRUE, 
                            check.names = FALSE, row.names = 1)

# Distance matrix
jaccard_mic_b <- hillpair(data = microbiota_bica, q = 1)
microbiota_jaccard_b <- as.dist(jaccard_mic_b$S)

# Absence/presence data frame - Jaccard distance
diet_jaccard_b <- vegdist(diet_bica, method = "jaccard")
diet_jaccard_b <- as.dist(diet_jaccard_b)


# Mantel test
set.seed(123)
mantel(diet_jaccard_b, microbiota_jaccard_b, method = "pearson",
       permutations = 999)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = diet_jaccard_b, ydis = microbiota_jaccard_b, method = "pearson",      permutations = 999) 
## 
## Mantel statistic r: 0.3353 
##       Significance: 0.006 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.135 0.182 0.219 0.284 
## Permutation: free
## Number of permutations: 999
# Make pcoas
dietPCoA_b <- as.data.frame(cmdscale(diet_jaccard_b))
plot(dietPCoA_b)

microbiotaPCoA_b <- as.data.frame(cmdscale(microbiota_jaccard_b))
plot(microbiotaPCoA_b)

# Procrustes analysis
procrust_b <- procrustes(X = dietPCoA_b, Y = microbiotaPCoA_b, scale=TRUE,
                         symmetric = TRUE)
pro_test_b <- protest(dietPCoA_b, microbiotaPCoA_b, permutations = 9999)
pro_test_b
## 
## Call:
## protest(X = dietPCoA_b, Y = microbiotaPCoA_b, permutations = 9999) 
## 
## Procrustes Sum of Squares (m12 squared):        0.7286 
## Correlation in a symmetric Procrustes rotation: 0.521 
## Significance:  0.0226 
## 
## Permutation: free
## Number of permutations: 9999
plot(pro_test_b)

eigen_b <- sqrt(procrust_b$svd$d)
percent_var_b <- signif(eigen_b/sum(eigen_b), 4)*100

beta_pro_b <- data.frame(procrust_b$X)
head(beta_pro_b)
trans_pro_b <- data.frame(procrust_b$Yrot)
head(trans_pro_b)
beta_pro_b$UserName <- rownames(beta_pro_b)
beta_pro_b$type <- "Diet"
beta_pro_b$seasons <- bica_metadata[,5]
beta_pro_b$species <- bica_metadata[,1]
head(beta_pro_b)
trans_pro_b$UsarName <- rownames(trans_pro_b)
trans_pro_b$type <- "Microbiota"
trans_pro_b$seasons <- bica_metadata[,5]
trans_pro_b$species <- bica_metadata[,1]
head(trans_pro_b)
colnames(trans_pro_b) <- colnames(beta_pro_b)
toplot_b <- data.frame(rbind(beta_pro_b, trans_pro_b))
pval_b <- signif(pro_test_b$signif, 1)

# graph diet and microbiota
procreustes_bica_new <- ggplot(toplot_b) +
  geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species, 
                                         shape = type)) + 
  scale_color_manual(values = c("#E69F00")) +
  theme_classic() +
  #scale_x_continuous(limits = c(-0.25, 0.4)) +
  #scale_y_continuous(limits = c(-0.2, 0.25)) +
  geom_line(aes(x = V1, y = V2, group = UserName), 
            col = "darkgrey", alpha = 0.6, 
            size = 0.2) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 10, colour = "black"),
        legend.position = "right",
        axis.text = element_text(size = 10, colour = "black"),
        axis.title = element_text(size = 12, colour = "black"),
        aspect.ratio = 1) +
  theme(legend.text = element_text(face = "italic")) +
  guides(color = guide_legend(ncol = 1)) +
  #annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_b), 
  #        size = 4) +
  xlab(paste0("PCoA 1 (",percent_var_b[1],"%)")) +
  ylab(paste0("PCoA 2 (",percent_var_b[2],"%)"))
print(procreustes_bica_new)

#ggsave("../figures/procreustes_bica_new.jpeg", width = 8.0, height = 4.0, dpi = 300)

Procrustes: Sceloporus grammicus

# Loading data

gram_metadata <- subset(metadata, Species_Scel == "Sceloporus grammicus")
diet_gram <- read.csv("../data/Proc_Diet_gram.csv", header = TRUE, check.names = FALSE, 
                      row.names = 1) %>% t()
microbiota_gram <- read.csv("../data/Proc_Mic_gram.csv", header = TRUE, 
                            check.names = FALSE, row.names = 1)

# Distance matrix
jaccard_mic_g <- hillpair(data = microbiota_gram, q = 1)
microbiota_jaccard_g <- as.dist(jaccard_mic_g$S)

# Absence/presence data frame - Jaccard distance
diet_jaccard_g <- vegdist(diet_gram, method = "jaccard")
diet_jaccard_g <- as.dist(diet_jaccard_g)

# Mantel test
set.seed(123)
mantel(diet_jaccard_g, microbiota_jaccard_g, method = "pearson",
       permutations = 999)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = diet_jaccard_g, ydis = microbiota_jaccard_g, method = "pearson",      permutations = 999) 
## 
## Mantel statistic r: 0.4457 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0905 0.1235 0.1568 0.1910 
## Permutation: free
## Number of permutations: 999
# Make pcoas
dietPCoA_g <- as.data.frame(cmdscale(diet_jaccard_g))
plot(dietPCoA_g)

microbiotaPCoA_g <- as.data.frame(cmdscale(microbiota_jaccard_g))
plot(microbiotaPCoA_g)

# Procrustes analysis
procrust_g <- procrustes(X = dietPCoA_g, Y = microbiotaPCoA_g, scale=TRUE,
                         symmetric = TRUE)
pro_test_g <- protest(dietPCoA_g, microbiotaPCoA_g, permutations = 9999)
pro_test_g
## 
## Call:
## protest(X = dietPCoA_g, Y = microbiotaPCoA_g, permutations = 9999) 
## 
## Procrustes Sum of Squares (m12 squared):        0.5908 
## Correlation in a symmetric Procrustes rotation: 0.6397 
## Significance:  1e-04 
## 
## Permutation: free
## Number of permutations: 9999
plot(pro_test_g)

eigen_g <- sqrt(procrust_g$svd$d)
percent_var_g <- signif(eigen_g/sum(eigen_g), 4)*100

beta_pro_g <- data.frame(procrust_g$X)
head(beta_pro_g)
trans_pro_g <- data.frame(procrust_g$Yrot)
head(trans_pro_g)
beta_pro_g$UserName <- rownames(beta_pro_g)
beta_pro_g$type <- "Diet"
beta_pro_g$seasons <- gram_metadata[,5]
beta_pro_g$species <- gram_metadata[,1]
head(beta_pro_g)
trans_pro_g$UsarName <- rownames(trans_pro_g)
trans_pro_g$type <- "Microbiota"
trans_pro_g$seasons <- gram_metadata[,5]
trans_pro_g$species <- gram_metadata[,1]
head(trans_pro_g)
colnames(trans_pro_g) <- colnames(beta_pro_g)
toplot_g <- data.frame(rbind(beta_pro_g, trans_pro_g))
pval_g <- signif(pro_test_g$signif, 1)

# Graph diet and microbiota
procreustes_gram_new <- ggplot(toplot_g) +
  geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species, 
                                         shape = type)) + 
  scale_color_manual(values = c("#56B4E9")) +
  theme_classic() +
  #scale_x_continuous(limits = c(-0.25, 0.4)) +
  #scale_y_continuous(limits = c(-0.2, 0.25)) +
  geom_line(aes(x = V1, y = V2, group = UserName), 
            col = "darkgrey", alpha = 0.6, 
            size = 0.2) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 10, colour = "black"),
        legend.position = "right",
        axis.text = element_text(size = 10, colour = "black"),
        axis.title = element_text(size = 12, colour = "black"),
        aspect.ratio = 1) +
  theme(legend.text = element_text(face = "italic")) +
  guides(color = guide_legend(ncol = 1)) +
  #annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_b), 
  #        size = 4) +
  xlab(paste0("PCoA 1 (",percent_var_g[1],"%)")) +
  ylab(paste0("PCoA 2 (",percent_var_g[2],"%)"))
print(procreustes_gram_new)

#ggsave("../figures/procreustes_gram_new.jpeg", width = 8.0, height = 4.0, dpi = 300)

Procrustes: Sceloporus spinosus

spi_metadata <- subset(metadata, Species_Scel == "Sceloporus spinosus")
diet_spi <- read.csv("../data/Proc_Diet_spi.csv", header = TRUE, check.names = FALSE, 
                     row.names = 1) %>% t()
microbiota_spi <- read.csv("../data/Proc_Mic_spi.csv", header = TRUE, 
                           check.names = FALSE, row.names = 1)

# Distance matrix
jaccard_mic_s <- hillpair(data = microbiota_spi, q = 1)
microbiota_jaccard_s <- as.dist(jaccard_mic_s$S)

# Absence/presence data frame - Jaccard distance
diet_jaccard_s <- vegdist(diet_spi, method = "jaccard")
diet_jaccard_s <- as.dist(diet_jaccard_s)


# Mantel test 
set.seed(123)
mantel(diet_jaccard_s, microbiota_jaccard_s, method = "pearson",
       permutations = 999)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = diet_jaccard_s, ydis = microbiota_jaccard_s, method = "pearson",      permutations = 999) 
## 
## Mantel statistic r: 0.1073 
##       Significance: 0.151 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.141 0.190 0.224 0.267 
## Permutation: free
## Number of permutations: 999
# Make pcoas
dietPCoA_s <- as.data.frame(cmdscale(diet_jaccard_s))
plot(dietPCoA_s)

microbiotaPCoA_s <- as.data.frame(cmdscale(microbiota_jaccard_s))
plot(microbiotaPCoA_s)

# Procrustes analysis
procrust_s <- procrustes(X = dietPCoA_s, Y = microbiotaPCoA_s, scale=TRUE,
                         symmetric = TRUE)
pro_test_s <- protest(dietPCoA_s, microbiotaPCoA_s, permutations = 9999)
pro_test_s
## 
## Call:
## protest(X = dietPCoA_s, Y = microbiotaPCoA_s, permutations = 9999) 
## 
## Procrustes Sum of Squares (m12 squared):        0.797 
## Correlation in a symmetric Procrustes rotation: 0.4505 
## Significance:  0.0789 
## 
## Permutation: free
## Number of permutations: 9999
plot(pro_test_s)

eigen_s <- sqrt(procrust_s$svd$d)
percent_var_s <- signif(eigen_s/sum(eigen_s), 4)*100

beta_pro_s <- data.frame(procrust_s$X)
head(beta_pro_s)
trans_pro_s <- data.frame(procrust_s$Yrot)
head(trans_pro_s)
beta_pro_s$UserName <- rownames(beta_pro_s)
beta_pro_s$type <- "Diet"
beta_pro_s$seasons <- spi_metadata[,5]
beta_pro_s$species <- spi_metadata[,1]
head(beta_pro_s)
trans_pro_s$UsarName <- rownames(trans_pro_s)
trans_pro_s$type <- "Microbiota"
trans_pro_s$seasons <- spi_metadata[,5]
trans_pro_s$species <- spi_metadata[,1]
head(trans_pro_s)
colnames(trans_pro_s) <- colnames(beta_pro_s)
toplot_s <- data.frame(rbind(beta_pro_s, trans_pro_s))
pval_s <- signif(pro_test_s$signif, 1)

# Graph diet and microbiota
procreustes_spi_new <- ggplot(toplot_s) +
  geom_point(size = 2, alpha = 0.75, aes(x = V1, y = V2, color = species, 
                                         shape = type)) + 
  scale_color_manual(values = c("#40A175")) +
  theme_classic() +
  #scale_x_continuous(limits = c(-0.25, 0.4)) +
  #scale_y_continuous(limits = c(-0.2, 0.25)) +
  geom_line(aes(x = V1, y = V2, group = UserName), 
            col = "darkgrey", alpha = 0.6, 
            size = 0.2) +
  theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 10, colour = "black"),
        legend.position = "right",
        axis.text = element_text(size = 10, colour = "black"),
        axis.title = element_text(size = 12, colour = "black"),
        aspect.ratio = 1) +
  theme(legend.text = element_text(face = "italic")) +
  guides(color = guide_legend(ncol = 1)) +
  #annotate("text", x = 0.07, y = -0.28, label = paste0("p-value=", pval_s), 
  #        size = 4) +
  xlab(paste0("PCoA 1 (",percent_var_s[1],"%)")) +
  ylab(paste0("PCoA 2 (",percent_var_s[2],"%)"))
print(procreustes_spi_new)

#ggsave("../figures/procreustes_spi_new.jpeg", width = 8.0, height = 4.0, dpi = 300)

Procustes and mantel analysis considering phylogenetic comosition of diet

tree.diet <- read.tree("../data/tree.nwk")
comm<-(diet)
IDs= colnames(comm)   # para obtener los nombres de las sp
tree.diet= ape::rtree(n=ncol(diet), tip.label = paste0(IDs))  # tener el árbol enraizado, con los nombres de los sitios y los nombres de las sp.


hill_phylo_parti_pairwise(comm,tree.diet, q = 0, show_warning = F) 
# Loading phylogenetic tree
phy_tree <- read_tree("../data/tree.nwk")

# Creating phyloseq object
SAM <- sample_data(metadata)
OTU <- otu_table(diet, taxa_are_rows=F)
phylo_physeq = phyloseq(OTU, SAM, phy_tree)

#Calculating phylogenetic distance
diet_dist_phy = distance(phylo_physeq, method = "uunifrac") %>% as.dist
set.seed(123)
# Order q=1
mantel(diet_dist_phy, microbiota_jaccard_1, method = "pearson",
       permutations = 999)

#Order q=2
mantel(diet_dist_phy, microbiota_jaccard_2, method = "pearson",
       permutations = 999)